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TECHNICAL  REPORT  SUMMARY 


The  principal  effort  during  the  period  covered  by  this  report  has  been 
directed  to  developing  more  thorough  understanding  of  the  proper  way  to  account 
for  the  effects  of  the  properties  of  the  real  Earth  on  the  amplitudes  and 
spectra  of  seismic  body  waves.  The  results  of  these  studies  will  contribute 
to  the  further  development  of  procedures  for  determining  the  energy  released 
at  the  source  and  the  mechanics  of  the  source  from  the  body-waves  observed 
at  large  distances. 

1.  A  full-wave  theory  applicable  to  the  synthesis  of  body-waves  at 
epicentral  distances  of  10°  to  40°,  for  propagation  through  an  upper  mantle  with 
either  first-order  discontinuities  or  intense  but  continuous  variations  in 
elastic  properties  has  been  developed.  The  theory  for  discontinuous  velocity 
changes  has  been  applied  to  the  synthesis  of  SH  waves  for  an  Earth  with 
discontinuities  in  density  and  velocity  at  depths  of  420  km  and  670  km.  Arrivals 
with  decreasing  amplitudes  at  distances  beyond  the  cusps:  predicted  by  geometric 
ray  theory  are  clearly  shown. 

One  way  of  modeling  a  region  of  rapidly  changing  properties  is  by  the  use 
of  a  large  number  of  thin  homogeneous  layers .  This  approach  requires  a 
lengthy  calculation  of  reflection-transmission  coefficients.  The  use  of  higher- 
order  Langer  approximations  of  the  radial  eigenfunctions  leads  to  a  more 
efficient  technique  for  treating  this  problem.  When  a  region  of  strong 
gradients  in  physical  properties  is  sufficiently  far  from  the  source  and 
receiver,  the  scattering  caused  by  the  region  results  in  waves  arriving  at 
points  far  from  the  receiver  or  at  times  removed  from  the  arrival  time  of  the 
transmitted  wave.  The  resulting  reduction  in  amplitude  of  P  and  SV  waves 
represents  a  strong  non-anelastic  attenuation. 
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2.  An  attempt  has  been  made  to  resolve  £  portion  of  the  frequency 

dependence  of  Q  ^  using  P  and  S  body-wave  spectra  from  shallow  and  deep 

earthquakes.  Various  theories  of  specific  dissipation  mechanisms  and  of  the 

anelastlc  attenuation  correction  of  body  waves  have  been  reviewed.  Under  the 

structure  of  those  theories,  a  model  of  Q  ( f )  is  proposed  in  which  Q  *  is 

constant  at  low  frequencies  but  decays  beyond  a  high  frequency  cutoff  which 

is  a  model  parameter.  The  model  is  tested  by  assuming  high  frequency  spectral 

-2  -3 

decay  slopes  of  both  to  and  .  P  and  S  waves  are  tested  independently, 

and  a  decay  in  the  absorption  spectrum  is  found  for  both.  However,  the 
required  modulation  is  not  the  same  for  the  two  wave  types,  suggesting  the 
existence  of  a  bulk-loss  mechanism  which  attenuates  a  narrow  range  of  P-wave 
frequencies.  The  depth  dependence  of  Q  *(f)  shows  that  the  bulk- loss 
mechanism  is  concentrated  almost  entirely  in  the  asthenosphere,  while  the 
shear  loss  mechanisms  are  distributed  in  a  complex  way  throughout  the  mantle. 

3.  Two  applications  of  elastodynamic  Green’s  functions  have  been  developed. 
First,  the  Green's  function  representation  of  a  wave  field  is  used  to  generate 
an  equivalent  elastic  field  that  analytically  reproduces  the  field  identical 

to  the  numerical  representation.  Then  this  analytical  representation  can  be 
used  easily  and  efficiently  to  predict  the  further  propagation  of  the  field  to 
long  distances.  The  second  application  uses  the  Green's  function  representation 
to  effect  a  transparent  grid  boundary  for  numerical  program  calculations. 


1.  BODY-WAVE  SYNTHESIS  AT  10°-40‘ 
V.F.  Cormier 


Upper  Mantle  with  First-Order  Discontinuities 


Theoretical  seismograms  enable  the  amplitude  and  waveform  of  body 
waves  to  be  incorporated  as  constraints  in  an  Inversion  scheme  for  an 
earth  model  or  the  source  time  function  of  an  earthquake  or  explosion. 

The  lower  mantle  has  long  been  known  from  travel  time  studies  to  be  nearly 
laterally  homogeneous  and  to  have  elastic  moduli  and  density  gradients 
that  vary  smoothly  and  slowly  over  large  ranges  of  depth.  Consequently 
in  the  distance  range  40-80°  simple  geometric  ray  theory  that  includes 
surface  reflections  combined  with  a  simple  source  description  is  sufficient 
to  synthesize  the  observed  waveform  of  body  waves  (Herrmann,  1975) .  At 
shorter  distances,  however,  geometric  ray  theory  must  be  abandoned  in 
favor  of  a  full  wave  theory  that  includes  non-ray  theoretical  effects  of 
waves  grazing  regions  in  the  upper  mantle  and  crust  having  discontinuous 
and/or  rapid  variations  in  velocity  and  density. 

Theoretical  seismograms  incorporating  such^  full  wave  theory  can  be 
generated  by  a  variety  of  methods,  but  in  general  all  start  with  a 
representation  of  the  form 


u(AQ,t)  -  ~  f  u,1'2  Re  ft  eiajJdpe_iWt  j 


(1.1) 


Eq.  1.1  is  appropriate  for  an  explosive  point  source  observed  at  the  earth's 
surface  at  distance  Aq  and  time  t  ,  where  f  is  a  path  in  the  complex 


ray  parameter  (p)  plane,  f 


is  a  product  of  reflection-transmission 
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coefficients  and  the  source-receiver  directivity  function,  and  J(p,A0) 
is  the  phase  delay  factor. 

The  phase  delay  factor  is  related  to  the  travel  time  and  distance 
Integrals  by  «• 


J(P»A0)  -  T(p)  -  pA(p)  +  pA0 


(1.2) 


More  specifically  the  phase  delay  factor  is  given  by 


J(p.Aq) 


rs 

/ 

rP 

/• 


[l/a2Cr)  -  p2/r2]1/2  dr 


[l/a2(r)  -  p2/r2] X^2  dr  +  pA0 


(1.3) 


where  rQ  is  the  radius  of  the  receiver,  rs  the  radius  at  the  source, 
rp  the  radius  at  which  the  integrand  vanishes,  and  the  integrand  is 
related  to  the  cosine  of  the  angle  of  incidence  by 


[l/a2(r)  -  p2/r2]1^2  -  cos  i/a(r) 


(1.4) 


Thus  the  radius  rp  is  the  turning  point  radius,  that  radius  where  the 
ray  bottoms  and  cos  i  =  0  . 

Fuchs  and  Muller's  (1971)  reflectivity  method  converts  the  integral 
over  the  ray  parameter  to  one  over  real  angles  of  incidence  in  which  the 
curvature  and  radial  heterogeneity  of  the  earth  is  modeled  by  a  stack  of 
homogeneous  flat  layers,  the  factor  f  representing  the  effect  of  the 
infinite  set  of  multiple  reflections  possible  in  the  layer  stack.  Some 
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methods  obtain  a  solution  in  the  time  domain  by  directly  operating  vith 
Eq.  1.1.  The  Cagniard  method  fora  sphere  (Gilbert  and  Helaberger ,  1972) 
achieves  this  by  operating  with  Eq.  1.1  in  terms  of  Laplace  transforms, 
in  which  -itu  -*■  s  ,  and  choosing  the  path  in  the  complex  ray  parameter 
plane  to  be  exactly  a  path  of  stationary  phase.  Chapman  (1976)  described 
how  the  solution  in  the  time  domain  can  also  be  obtained  by  the  first 
motion  approximation,  which  evaluates  the  double  integral  in  Eq.  1.1 
by  the  equal  phase  method. 

Richards  (1973)  enumerated  the  advantages  that  accrue  to  the  more 
general  procedure  of  solving  the  inner  integral  in  the  frequency  domain 
by  a  numerical  integration  along  a  path  T  in  the  complex  ray  parameter 
plane  and  then  obtaining  the  time  domain  solution  by  inverse  Fourier 
transformation.  In  review,  these  advantages  are  (i)  that,  unlike  the 
Cagniard  method,  the  path  T  may  remain  fixed  for  a  series  of  distances, 
the  only  constraints  being  that  T  be  sufficiently  near  all  ray  theory 
saddles  and  end  in  regions  in  which  the  integrand  is  exponentially  small; 

(ii)  that  the  method  can  be  extended  for  ray  paths  having  a  turning  point; 

(iii)  that  the  path  T  automatically  includes  the  non-ray  theoretical 
effects  of  diffraction  given  by  the  residue  contribution  of  ray  parameter 
poles  in  the  reflection-transmission  coefficients;  and  (iv)  that  effects 
of  attenuation  are  easily  incorporated  by  allowing  the  velocity  profile 
to  be  complex. 

Incorporation  of  Attenuation 

Since  the  solution  formulae  (Eqs.  1.1-1. 4)  are  all  analytical  in  velocity, 
the  solution  for  an  attenuating  earth  is  given  by  analytic  continuation 
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to  a  complex  velocity  profile  defined  in  terms  of  a  Q  and  a  real 

velocity  model  (Cormier  and  Richards,  1976).  The  observation  that  the 

Q  of  the  earth  is  observed  to  be  nearly  frequency  independent  over  the 

-2  3 

broad  band  of  seismic  frequencies  (10  -  10  Hz)  coupled  with  the 

hypothesis  that  the  attenuation  mechanism  be  linear  leads  to  the  relation 
for  the  complex  velocity  profile 


v(r,w2) 


VR^OJi)  1  + 


hQqOO 


&n 


(“i)‘  2Vr>  ] 


(1.5) 


(Liu  et  al. ,  1976),  where  vR(r,u)^)  is  the  real  velocity  profile  determined 
at  the  reference  frequency  and  QQ(r)  is  the  Q  profile.  Many 

other  representations  for  v  can  be  derived  by  allowing  for  different 
distributions  of  discrete  relaxation  mechanisms  with  unequal  strengths 
over  the  same  frequency  band. 

Under  several  conditions  the  dispersive  term  in  Eq.  1.5  can  be 
neglected,  resulting  in  the  relation 


v(r)  -  vR(r>  [  1  ~  2Q^T  ]  (1'6) 

The  conditions  under  which  dispersion  can  be  neglected  are  that 

1.  calculations  incorporating  v(r)  be  over  a  sufficiently  narrow 
frequency  band; 

II.  Q0(r)  be  sufficiently  large;  and 

the  real  velocity  profile  vR(r)  be  determined  from  data  in  a 
frequency  band  coincident  with  that  of  the  calculation. 


III. 


When  either  condition  I  or  II  is  not  met,  the  dispersive  effect  of 
attenuation  will  be  manifest  in  a  change  of  the  shape  and  rise  time  of 
a  propagating  pulse  (Futterman,  1962).  Although  conditions  I  and  II 
affecting  pulse  shape  may  be  satisfied  by  the  pass  band  of  most  seismograph 
systems  combined  with  reasonable  models  of  the  anelasticity  of  the  upper 
mantle,  condition  III,  which  can  bias  travel  times,  may  not.  Thus 
knowledge  of  the  deviation  from  "average"  anelastic  properties  along  a 
ray  path  may  assist  in  applying  a  travel  time  correction  when  the  arr 
time  of  that  ray  is  used  in  a  location  determination. 

A  time  domain  study  of  anelasticity' s  effects  on  amplitudes,  tra 
times,  and  rise  times  of  body  waves  may  clarify  the  nature  of  the  frequency 
dependent  anelasticity  inferred  from  body  wave  spectra  as  reported  by 
Lundquist  (1976). 

Preliminary  Results 

The  method  described  by  Cormier  and  Richards  (1977)  for  synthesizing 
the  seismogram  of  a  body  wave  interacting  with  a  discontinuous  velocity 
increase  has  been  extended  to  synthesize  seismograms  interacting  with 
two  or  more  such  increases  with  any  arbitrarily  close  spacing  in  depth. 

A  test  of  the  method  was  begun  by  synthesizing  the  waveform  of  a  SH 
wave  having  a  delta  source  time  function  interacting  with  discontinuous 
velocity  and  density  increases  at  420  and  670  km  depth.  Figure  1.1  is  a 
reduced  travel  time  curve  for  SH  waves  in  this  earth  model.  Initial 
calculations  Incorporated  a  complex  velocity  profile  of  the  type  given 
by  Eq.  (5)  but  with  QQ(r)  taken  to  be  1  10^  .  Later  calculations  will 


test  the  effect  of  anelasticity  on  waveform  and  apparent  travel  times 


Fig-  1.1  Reduced  travel  time  curve  for  SH  waves  in  the  upper  mantle  of 
earth  model  1066B.  Weaker  arrivals  (x's)  at  distances  less 
that  E  and  C  are  partial  reflections  from  discontinuous  velocity 
increases  at  420  and  670  km.  At  distances  greater  than  20°  the 
x's  represent  the  interference  head  waves  along  these  discontinuities. 
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xtr 


with  reasonable  Q  models  of  the  upper  mantle. 

Figure  1:2  shows  the  synthesized  SH  waveforms  in  the  distance  range 
30-34°.  Mote  that  the  abrupt  cutoff  of  a  ray  arrival  at  cusp  D  at 
28°  predicted  by  geometrical  ray  theory  (valid  at  infinite  frequency) 
becomes  (at  finite  frequency)  a  gradual  decay  in  amplitude  at  distances 
greater  than  D  .  Mote  also  that  the  interference  head  waves  associated 
with  the  420  and  670  km  discontinuities  persist  for  long  distances  at 
which  the  simply  transmitted  ray  bottoms  below  both  discontinuities. 

The  anqtlitude  of  such  head  waves  has  been  shown  to  be  strongly  sensitive 
to  the  anelastic  properties  along  the  underside  of  the  velocity 
discontinuities  (Hill,  1971). 

Interference  of  the  travel  time  branches  results  in  a  maximum  peak— 
to-peak  amplitude  of  the  SH  body  wave  near  20“  (Helmberger  and  Engen, 
1974).  The  amplitude  growth  near  20°  will  be  diagnostic  of  both  the 
elastic  and  anelastic  properties  of  the  earth  in  the  depth  range  300-700 
km.  Only  a  synthesis  incorporating  both  properties  as  depth  dependent 
parameters  can  separate  their  competing  effects  (Kennett,  1975). 

Synthetics  for  long  period  SH  body  waves  will  be  completed  for  the 
distance  range  12°-34°  for  an  earth  model  with  sharp  discontinuities  in 
the  upper  mantle  and  for  several  attenuation  mechanism  models  in  the  depth 
range  100-700  km.  It  is  hoped  that  by  comparison  with  observed  SH 
waveforms  in  this  distance  range  some  new  constraints  may  be  obtained 
regarding  elastic  and  anelastic  structure  in  the  upper  mantle. 
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Fig-  1.2  Response  of  earth  model  1066B  to  a  delta  source  time  function 


in  SH  displacement  (surface  focus)  convolved  with  the  instrument  response 


of  a  15-100  WWSSN  seismograph.  (1)  SH  wave  transmitted  through 


670  km  discontinuity  plus  interference  head  wave  along  that 


discontinuity;  (2)  total  reflection  plus  diffraction  from  420 


km  discontinuity ;  (3)  interference  head  wave  associated  with 


420  km  discontinuity 
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Upper  mantle  with  regions  of  Intense  but  continuous  variations 
in  elastic  properties 


Introduction 

By  proper  choice  of  the  radial  eigenfunctions  that  solve  the 
potential  wave  equation  in  an  inhomogeneous  medium,  the  number  of  layers 
required  to  describe  the  medium  is  minimized  and  thereby  also  the  number 
of  reflection-transmission  coefficients  required  to  be  evaluated  in  the 
factor  f  in  Eq.  1.1.  By  allowing  the  layers  to  be  inhomogeneous,  in 
which  the  elastic  properties  vary  slowly,  only  the  coefficients  associated 
with  first-order  discontinuities  need  to  be  evaluated.  Langer's  uniformly 
asymptotic  approximation  to  the  radial  eigenfunctions  can  be  used  to 
evaluate  coefficients  that  are  valid  for  ray  parameters  corresponding  to 
rays  both  near  grazing  and  steep  incidence  to  a  discontinuity  (Richards, 
1976). 

Suppose  now,  however,  that  the  medium  contained  no  first  order 
discontinuities  but  only  thin  regions  in  which  elastic  properties  varied 
rapidly,  bounded  by  at  most  second-order  discontinuities.  Figure  1.3 
compares  two  possible  velocity  distributions  in  the  upper  mantle 
assuming  either  first  order  discontinuities  or  continuous  variation  in 
elastic  properties.  A  representation  of  displacement  such  as  Eq.  1.1 
can  be  obtained  for  a  continuous  medium  in  either  the  Cagniard  or 
reflectivity  methods  by  modeling  the  regions  of  rapid  variation  by  a 
sequence  of  many  thin  homogeneous  layers. 
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Second  Order  Wave  Equations 

To  avoid  a  lengthy  calculation  for  the  factor  f  to  Include  the 
many  possible  multiple  reflections  In  such  a  stack  of  layers,  one  can  seek 
a  higher— order  approximation  to  the  radial  eigenfunction.  Such  an 
approximation  must  be  accurate  at  the  lowest  frequencies  required  to 
synthesize  long  period  body  waves  and  can  be  obtained  from  the  second- 
order  Helmholtz  equations  satisfied  by  P,  SV,  and  SH  potentials.  These 
potentials  are  given  by  Richards  (1974)  as 


y2p  +  fep  +  v 


V2V  + 

y 


_  fHi  +  o/M\ 

1/2.2  3r  (w2  ) 


(1.7a) 


p  0) 
K 


V  +  ESVV  "  1/2 

p  in 


1  _  /„  _i_  ^  ifr\  +  0l2l)a 

2  2  l  0  ’  r  Bin  6  2<f>  »  r  30  )  2}  V 

0)  \  T  f  U)  / 


v2h  +  H  +  etH 


(1.7c) 


for  P,  SV,  and  SH  waves,  respectively.  (1  ,  p  are  elastic  constants, 

p  density,  et  ,  egv  ,  ep  ,  functions  of  the  density,  elastic 

constants,  and  their  first  and  second  order  radial  derivatives,  and 

U  =  (u  ,  Uo  ,  Uj.)  is  the  vector  displacement.) 
r  o  <J> 

A  separation  of  variables  allows  the  angular  dependence  of  P  ,  y  * 
and  H  to  be  satisfied  by  vector  surface  harmonics  that  can  be  factored 
out  of  Eqs.  1.7a-1.7c,  resulting  in  the  scalar  radial  equations 


d2x 


dr 


f  +  u>2Q2  X  +  ep  X 


d2 (X/r)  B9Y 

2  “  2 
.  dr  r  3r 


.7b) 


(1.8a) 
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~  +  to2Q2  Z  +  et  Z  *  0 


(1.8c) 


where  Q 


1  "  (a2" 

2  '  (?' 


u  2  2 

B  =  -  0)  p 


P2/r2N|  1/2 


p2/r2 )  1/2 


for  a  P  velocity  function  a  ,  S  velocity  function  8  ,  and  ray 
parameter  p  .  The  functions  X  ,  Y  ,  and  Z  divided  by  r  are  the  radial 
eigenfunctions  of  the  full  potentials  P  ,  V  ,  and  H  ,  respectively. 


First-Order  Solutions 

Within  a  constant  of  proportionality  to  be  determined  by  a 
normalization  criterion,  Langer's  (1951)  uniform  approximation  to  the 
solutions  X  ,  Y  ,  and  Z  can  be  expressed  at  the  lowest  order  in  frequency 
as  (Richards,  1976) 


(1) 

1/2  „(2) 


X  -  (  Ki/Qi  ) 11  H<2'  (a%) 


(1.9a) 


Y  -  (c2/q2  Jl 1/2  h{/3  («e2)/(up 


(1.9b) 


z  =  (  S2/Q2  )  1/2  *ul  <*2> 


(1.9c) 


where  is  a  spherical  Hankel  function  of  order  1/3  and  kind  (1)  for 
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an  upgolng  or  (2)  for  a  down  going  wave  and  2  are  defined  by 


*  /  V 


(1.10) 


In  order  Chat  the  coupled  equations  for  P  and  V  (X  and  Y)  be  dimensionally 
correct,  the  factor  ujp  (constant  in  r)  must  appear  as  shown  in  the 
solution  for  Y  . 

To  first  order  in  frequency  the  P  and  SV  potentials  are  decoupled 
and  the  properties  of  the  first  and  higher  order  radial  derivatives  of 
elastic  moduli  and  density  given  by  the  factors  e  ,  ec  ,  and  e„  do 

P  bV  1 

not  appear  in  the  solution.  The  solutions  given  above,  however,  do 
correctly  account  to  first  order  for  the  effects  of  radial  inhomogeneity 
near  the  turning  point  of  a  ray.  These  solutions  are  valid  for  body 
waves  of  frequency  0.01  Hz  or  more  except  in  regions  of  severe  velocity 
and/or  density  gradient  in  the  upper  mantle  low  velocity  and  transition 


Higher  Order  Solutions 

Higher  order  uniform  approximations  may  be  determined  from  a 
perturbation  solution  for  the  functions  A(r,u>)  and  B(r,a>)  in 


V  °  A(r,co)u  +  B(r,u))u' 


(1.11) 


where  V  symbolizes  X  ,  Y  ,  or  Z  and  u  symbolizes  a  first  order  solution 
given  by  Eq.  1.9a,  b,  or  c.  Earlier  studies  (Cherry,  1950;  Olver,  1974; 
Chapman,  1976)  work  with  the  Langer  transformation  of  Eqs.  1.8a-c  which 


for  a  particular  radial  function  V  is  of  the  form 


0  (1.12) 

,  *  -  6(r)  -  -  -  (3/2C2/3) 

e 

A  perturbation  solution  seeks  the  evaluation  of  the  functions  A(r9u>)  and 
B(r,oj)  in 

V  =  A(z,uj)  bi(oj2/3z)  +  B(z,<d)  bi'(u)2/3z)  (1.13) 

where  bi  denotes  an  Airy  function  of  the  first  or  second  kind  and  (') 
a  derivative  with  respect  to  z  .  It  should  be  noted  that  the  solution 
given  by  Eq.  1.13  is  equivalent  to  that  given  by  Eq.  l.n  after 
back-transforming  V  to  V  and  expressing  bi  in  terms  of  Hankel  functions 
of  order  1/3  . 

The  evaluation  of  the  functions  A  and  B  involve  integrals  over  the 
variable  z  of  integrands  containing  the  function  a>g(z,h))  .  Numerical 
Integration  over  z  in  seismic  applications  requires  that  simple  functions  of 
radius  in  tug(z ,w)  such  as  e  be  determined  as  complicated  functions 
of  z  .  If  the  variable  of  integration  were  changed  from  z  to  r  ,  functions 
such  as  e  could  be  easily  evaluated,  but  the  functions  0*  ,  0"  ,  and 
0'"  resulting  from  the  Langer  transformation  would  be  time 
consuming  to  evaluate  as  functions  of  r  .  Thus  the  alternative  solution 
form  given  by  Eq.  1.11  will  be  sought  because  the  integrals  needed  to 
determine  the  coefficients  of  A(r,u>)  and  B(r,u>)  in  power  series  of 


_  [o)2z  _  g(z  ,to)  ]  V 
dz 


where  ug(<.u)  -  +  f  ^4 


1  d 

2 


and  V 


0 


’1/2v  . 


e1. 
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u 


of  frequency  w  are  Integrals  over  radius  of  simple  functions  of  radius. 
In  seeking  such  a  solution,  let 


where 


X  -  A(r,uj)x  +  B(r,w)x’ 


Y  =  X(r,oi)y  +  B(r,0))y' 


Z  -  C(r,u>)y  +  D(r,oj)y' 


x  -  (Sl/Ql)1/2  H^O*)^) 


y  =  (C2/Q2) 1/2  H^/3(a)52)/a>p 


y  =  <C2/Q2>1/2  Hi/3(^2) 


A(r,w)  +  u“2nAn(r) 
n-0 


B(r,w) 


(1.14a) 

(1.14b) 

(1.14c) 

(1.15a) 

(1.15b) 

(1.15c) 

(1.16a) 

(1.16b) 


n«0 


and  similarly  for  A  ,  B  ,  C  ,  and  D  .  For  example,  in  ascending  negative 
powers  of  frequency  the  first  three  terms  in  the  solution  for  Z  would  be 

1/2  i  **2/3  (^2) 

2  *  Vg2/Q2)  Hi/3<^2>  +  ~ —  '  ~m  —  -  2 


cx  (C2/Q2)1/2  +  d1(C2/Q2)“1/2 


Hi/3^2> 


w 


(1.17) 
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To  determine  the  coefficients  A  ,  B  ,  etc. »  the  solution  forms  of 

n  n 

Eqs.  I.l4a-c  are  substituted  in  the  radial  Eqs.  1.9a-c  and  terms  in  equal 
powers  of  frequency  are  equated. 

The  results  for  the  zeroth  and  first  higher  order  term  for  P  ,  SV  , 
and  SH  are  as  follows: 


X 


ojpY 


Z 


a  .  1/2  j  ,  -  .  LjIP^i/Ql)1/  H1/3(U)€1) 

VVV  H1/3(“V  +  2o)  +  * 

—  1/2  1  **1*SV^2^2^  ^2/3  (^2) 

V  W  hJ1/3^2>  +  J-Jy-— i ^  + 


(1.18a) 


.  (1.18b) 


hj(?2/Q2)1/2  h{/3(0)52)  + 


L1ISH^2/Q2)  /  H2/3^2* 


.  (1.18c) 


where  and  are  normalization  constants  for  up-  (j=l)  or  down- 

( j=2)  going  waves.  The  functions  Ip  ,  Igv  ,  Igy  are  defined  by 


(1.19a) 


(1.19b) 


(1.19c) 


and  the  coupling  coefficients 
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and 

The  coupling  coefficients  as  they  appear  in  the  second  order  term  for  P  and 
SV  waves  represent  singly  converted  waves  that  ususally  arrive  at  times 
and  distances  removed  from  those  of  a  wave  directly  transmitted  through  a 
region  of  anomalously  large  gradient  when  the  region  is  much  deeper  than 
either  the  source  or  receiver. 

To  accurately  calculate  the  third  order  term  for  the  P  and  SV 

2 

potential  solutions,  the  corrections  0(X/w)  and  0(Y/w  )  noted  in 

Eqs.  7a-b  must  be  found.  The  results  of  calculating  the  next  higher 

order  term  for  SH  indicate  that  this  term  is  of  order  Ig^y/w  .  Thus 

we  may  predict  that  the  third  order  terms  for  P  and  SV  are  of  order 
2  2  2  2 

I pX/w  and  I^y/w  .  The  third  order  P  and  SV  terms  will  also  involve 

coupling  coefficients  appearing  in  double  integrals  of  the  form 


\  2  li2/3t“Sl) 

/  r  hJ.JwL) 


Physically  these  double  integrals  involving  coupling  coefficients  represent 
waves  that  travel  essentially  as  converted  waves  in  regions  of  intense 
gradient,  converting  back  to  the  original  wave  type  after  leaving 
the  region  on  their  way  towards  the  receiver. 
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Displacement 

The  higher  order  terms  in  the  Langer  approximation 
may  be  applied  to  determining  the  reflection-transmission  coefficients 
of  elastic  waves  incident  on  a  radial  discontinuity  of  first,  second, 
or  third  order  in  elastic  moduli  or  density.  Thus  fewer  but  more 
complicated  factors  need  to  be  evaluated  in  the  determination  of  the 
factor  f  in  Eq.  l.i. 

Let  us  now  instead  consider  the  problem  of  wave  propagation  in  a 
hypothetical  medium  in  which  the  elastic  moduli  and  density  are  everywhere 
analytic  in  depth  but  exhibit  intense  gradients  in  certain  depth  ranges. 
Although  elastic  moduli  and  density  and  their  first  and  higher  order 
derivatives  may  rapidly  change  in  certain  zones  in  the  earth,  let  us 
assume  that  they  nevertheless  change  continuously.  The  advantages  given 
by  this  assumption  are  (1)  that  the  representation  for  potential  or 
displacement  at  the  surface  as  an  integral  over  ray  parameter  does  not 
require  reflection-transmission  coefficients  to  appear  in  the  integrand, 
and  (2)  that  the  WKBJ  approximation  to  the  radial  eigenfunctions  may  be 
substituted  because  these  functions  need  only  to  be  evaluated  at  the 
source  and  receiver  radius  far  from  a  turning  point. 

Using  the  equations  relating  displacement  U  and  potential 


U  =*  u1/2(o,0,-  )  for  SH  waves  and 


U20) 


li 


7  grad 
1 


(~I/2  -T-)  +f  curl  Curl(-I72  7  •  0  •  °) 


K-i  K«  n 

-ij  (U _  ,  0  ,  0)  +  y- 
pui  pu> 


(1.  21) 
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for  P  and  SV  waves,  a  representation  for  the  Fourier-transformed 
displacement  as  an  integral  over  ray  parameter  may  be  obtained  as  outlined 
by  Richards  (1973).  Including  now  the  first  higher  order  term  in  the 
potentials  results  in  a  displacement  representation  for  a  source  at  the 
receiver  radius  and  a  ray  departing  downwards  from  the  source: 

r  .  irc/2_ 

Ur(r,A)  -  o)1/2  Lp  J  P1/2  ^1  + - —  £^e*U>J(r,p>  dp  (1.22a) 

r 

Ue(r.A)  -  J-'1  Lsyft12  (l  *—=~Si)  ^(r-P>  OP  (l-22b) 

r 

r  Ur/2 

U,(r.A)  -  J'2  LS„  / ^  (l  ♦  dp  (1.22c) 

r 


when  only  the  radial  displacement  Ur  of  a  P  wave  and  horizontal 

VW  1 

displacement  Ug  of  a  SV  wave  is  calculated.  The  constants  Lp  ,  Lgy  , 

Lou  are  proportional  to  the  source  normalization  factor  and  seismic 

moment.  The  phase  factor  J  is  calculated  in  the  P  velocity  profile  and 

J  in  the  S  velocity  profile.  Note  that  the  tt/2  phase  shift  of  the 

higher  order  term  in  the  integrand  may  be  interpreted  as  the  phase  shift 

expected  of  reversed  travel  time  branches  generated  by  regions  of  intense 

gradient. 

Summary  and  Conclusions 

A  higher  order  Langer  approximation  has  been  developed  for  the  radial 
eigenfunction  of  the  second  order  wave  equations  satisfied  by  P  and  S 
potentials.  The  higher  order  approximation  allows  a  more  exact  solution 
to  the  elastic  wave  equations  in  regions  of  strong  gradient  of  velocity 
and  density,  remains  valid  for  a  turning  point  in  such  regions,  and 


~9 
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minimizes  the  number  of  radially  inhomogeneous  layers  required  to  describe 
an  earth  model.  Coefficients  in  the  perturbation  series  solution  are 
expressed  as  integrals  of  functions  of  velocity,  density,  and  their 
radial  derivatives  over  the  depth  coordinate  rather  than  over  the  Langer 
transformed  coordinate.  The  coefficients  can  be  inexpensively  evaluated 
to  third  order  in  frequency  for  SH  and  to  second  order  for  P  and  SV  waves. 
The  second  order  term  in  the  approximation  for  P  and  SV  waves  includes 
the  effect  of  single  P-SV  scattering  by  regions  of  strong  gradient.  When 
such  regions  are  stuff iciently  far  from  the  source  and  receiver,  these 
scattered  waves  arrive  either  at  distances  far  from  the  receiver  or  at 
times  removed  from  the  arrival  time  of  the  transmitted  wave.  Thus  they 
may  represent  an  important  source  of  non- intrinsic  attenuation  for  low 
frequency  P  and  SV  waves  in  the  upper  mantle. 

Test  calculations  are  now  in  progress  using  the  displacement 
representation  given  by  Eqs.  1.22a-c  in  simple  models  of  upper  mantle 
transition  zones.  Comparisons  of  seismograms  calculated  using  this 
representation  for  waves  interacting  with  a  thin  zone  of  rapidly  increasing 
velocity  and  density  will  be  made  with  seismograms  interacting  with  a 
discontinuous  velocity  increase  at  the  same  depth.  These  comparisons  and 
calculations  will  be  used  to  clarify  how,  if  at  all,  the  representation 
given  by  Eqs.  1.22a-c  may  acoount  for  narrow  as  well  as  wide  angle 


reflections  from  transition  zones. 
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2.  THE  FREQUENCY  DEPENDENCE  OF  Q 
Gary  M.  Lundquist 

Introduction. 

Energy  absorption  in  the  Earth  is  generally  measured  in  terms  of  the 
seismic  quality  factor,  Q  ^  ,  which  is  defined  as  the  fraction  of  energy 
dissipated  per  wavelength:  Q  *  *  —  ~  .  Since  Q  *  is  defined  on  a  per 
cycle  basis,  then  a  constant  Q  ^  implies  increasing  attenuation  with 
frequency  for  a  given  length  of  path.  In  particular,  anelastic  attenuation 
is  low  at  low  frequencies  (f  <  0.1  Hz)  ,  amounting  to  a  few  percent  at  most 
over  teleselsmlc  path  lengths.  Mapping  Q  ^  as  a  function  of  frequency  at 

those  wavelengths  requires  a  method  with  high  resolution.  On  the  other  hand, 
anelastic  attenuation  for  f  >  1  Hz  may  reduce  wave  amplitudes  by  several 

orders  of  magnitude,  so  that  a  low  resolution  technique  may  be  used  to 
observe  Q  ^(f) 

The  functional  dependency  of  Q  ^  on  frequency  has  been  the  subject  of 
many  recent  papers  with  two  basically  different  thrusts.  The  first  group 
presents  observational  evidence  on  Q  1  and  concludes  that  Q_1  is  independent 
of  frequency  over  the  period  range  from  about  1  hour  to  1  sec.  This  literature 
has  been  concisely  reviewed  by  Anderson  and  Hart  (1977),  and  includes  data 
from  free  oscillations,  surface  waves,  and  body  waves  with  emphasis  on 
multiple  ScS  paths.  Most  of  the  accurate  determinations  of  Q  *  structure 
or  average  Q  ^  over  teleseismic  paths  are  done  for  waves  with  periods  greater 
than  5  sec,  which  is  also  a  period  range  of  small  attenuation  and  therefore 
low  resolution. 

A  second  class  of  papers  has  examined  the  physical  dissipation  mechanisms 
which  could  be  responsible  for  anelastic  attenuation.  These  papers  were 


thoroughly  reviewed  by  Jackson  and  Anderson  (1970),  who  find  a  marked  frequency 
dependence  for  Q  ^  associated  with  each  individual  mechanism.  One  of  the  more 
likely  mechanisms  for  anelastic  attenuation  is  grain-boundary  relaxation 
(Solomon,  1971;  Anderson  and  Hart,  1977),  which  may  be  modelled  as  a  relaxation 
process  in  a  standard  linear  elastic  solid  (Mason,  1958).  This  model  predicts 
Q  ^  to  behave  as  to  for  low  frequencies  and  as  to  *  for  high  frequencies 
relative  to  the  "relaxation  frequency"  at  which  absorption  is  a  maximum.  Of 
course,  grain  sizes,  existence  and  viscosity  of  melt  phases,  and  varying  activation 
energies  will  lead  to  a  range  of  relaxation  times  and  therefore,  a  page  of 

absorption  peaks,  yielding  a  weaker  frequency  dependence  than  would  be  expected 
from  a  single  mechanism. 

Liu  et  al.  (1976)  chose  a  specific  distribution  of  relaxation  times  in 
a  standard  linear  elastic  solid  which  gave  them  a  frequency  independent 

Q  1  over  a  frequency  band  wider  than  the  band  of  observation.  Their  model 
may  be  visualized  as  an  absorption  spectrum  with  a  finite  bandwidth  and 
to  and  to  1  decays  at  lower  and  higher  frequencies  respectively.  The  existence 
of  a  low  frequency  cutoff  is  required  for  such  a  model,  or  the  real  part  of 
the  index  of  refraction  would  be  unbounded  at  zero  frequency  (Futterman,  1962), 
but  the  position  of  the  cutoff  is  not  determined  either  by  Futterman's  theory 
or  by  Liu  et  al.  In  the  absence  of  evidence  to  the  contrary,  the  low 
frequency  cutoff  was  placed  beyond  the  range  of  observation  at  10~^  Hz. 

Again,  no  presently  available  technique  has  sufficient  resolution  to  resolve 
the  existence  of  this  cutoff,  much  less  its  position  in  frequency. 

Anderson  and  Hart  (1977)  estimate  the  total  bandwidth  of  the  absorption 
spectrum  by  using  other  theories  to  estimate  the  total  velocity  dispersion, 
which  is  a  required  consequence  of  any  constant  Q~^  model  (Newlands,  1954). 

That  is,  if  the  difference  in  phase  velocity  at  the  low  and  high  frequency 
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cutoffs  can  be  independently  estimated,  then  the  bandwidth  is  predicted  by 
the  model  of  Liu  et  al.  Anderson  and  Hart  estimate  a  15  to  20Z  change  in 
velocity,  which  corresponds  under  certain  assumptions,  to  a  bandwidth  of 
3-6  decades.  If  the  low  frequency  cutoff  is  at  10  ^  Hz,  then  five  decades 
gives  a  high  frequency  cutoff  at  10  Hz,  which  is  used  by  Liu  et  al.  (1976). 

A  high  frequency  cutoff  is  required  by  causality  (Lamb,  1962)  but  is  not 
uniquely  determined  either  by  causality  or  by  Anderson  and  Hart.  The  object 
of  this  report  is  to  present  evidence  for  a  high  frequency  cutoff  within  the 
passband  of  WWSSN  seismographs. 

The  present  study  of  anelastic  attenuation  arose  out  of  a  study  of 
body-wave  seismic  source  spectra.  The  standard  correction  for  the  Earth's 
attenuation  filter  is  (Ben  Menshem  et  al.,  1966) 

G(x)exp(irf  T/Qg)  -  G(x)exp(irft*) 

where  G(x)  is  the  frequency  independent  geometric  attenuation.  Evaluation  of 
the  travel  time,  T  ,  and  the  average  or  "effective"  Q  ^  ,  requires  both 

elastic  and  anelastic  Earth  models.  Seismic  velocities  are  well  constrained, 
but  Q  1  models  are  still  in  a  state  of  flux. 

Of  the  models  available,  two  California  Institute  of  Technology  Q_1 
models  were  tested  on  spectra  from  several  central  Asian  crustal  earthquakes. 

CIT  208  was  derived  from  short-period  body  waves  (Julian,  personal  communication), 
and  for  shallow  events  could  be  represented  by  a  distance  independent 

■k 

tg  ~0.42  .  CIT  11  CS2  -  QM  ,  on  the  other  hand,  was  determined  from 

surface  wave  data  (Anderson,  1967)  and  could  be  represented  by  t*  -  1.0 
Of  course  there  is  no  explicit  frequency  dependence  in  either  model. 


The  discrepancy  in  t  could  be  interpreted  either  as  an  error  in  one 

*  -1 

model  or  as  evidence  that  t  (or  Q  )  is  frequency  dependent.  The  first 

model  gave  intuitively  acceptable  corrections  to  the  body  wave  spectra;  while 

the  second  model  so  overcorrected  the  high  frequencies  that  spectral  content 

at  the  peak  of  the  short-period  WWSSN  seismometer  response  would  have  to  be 

discarded  as  noise  (see  Figure  2.3).  Though  the  first  model  gave  usable 

spectra,  the  published  literature  (Carpenter  and  Flynn,  1965;  Langsten  and 

Helmberger,  1975)  suggested  the  second  model  was  more  correct. 

The  conflict  was  resolved  by  assuming  that  Q  decayed  toward  high 

frequency  according  to  a  standard-linear-elastic-solid  model.  By  assuming 

spectral  shape,  and  in  particular  by  requiring  the  high  frequencies  to  decay 
-2  -3 

as  Cl)  or  a)  ,  the  position  of  the  high  frequency  decay  could  be  observed 
for  the  raypaths  tested.  At  this  time  it  is  not  possible  to  determine 
whether  this  decay  is  the  high  frequency  cutoff  postulated  by  Liu  et  al. 
(1976),  or  simply  a  window  in  an  otherwise  constant  absorption  spectrum. 
Theoretical  Background. 

The  theory  will  be  briefly  reviewed  to  emphasize  the  necessary  model 
parameters.  The  energy  of  an  harmonic  wave  may  be  written  in  terms  of  its 
initial  value,  Eq(u))  >  and  a  spreading  factor  G(x)  ,  as 

E(w)  =  G  E  (u>)  2i(Kx  "  Wt)  (2.1) 

o  e 

where  K  *  +  iic' '  is  a  complex  wave  number.  Q  ^  has  already  been  defined 

as  the  fraction  of  energy  dissipated  per  wave  length ,  X  , 

n-i  _  _1_  E(nX)  -  E((n  +  1)A) 

W  2tt  E(nX) 


(2.2) 
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If  energy  loss  Is  small,  then  to  first  order  the  energy  carried  is  given 
by  2.1  and 


-2tc' 

e 


For  small  kc"  ,  the  series  expansion  of  the  exponential  gives 


(2.3) 


That  is,  Q~*  may  be  defined  in  terms  of  the  ratio  of  the  imaginary  and  real 
parts  of  the  wave  number. 

The  distinction  between  P-wave  and  S-wave  attenuation  is  obtained  by 
writing  the  respective  wave  numbers  in  terms  of  complex  velocities. 


o  2  2 

2  lo 

*6  "  g2  "  p/p 


where  K  and  U  are  complex  bulk  and  shear  moduli  respectively.  For  the 
case  of  small  attenuation,  the  ratio  of  imaginary  to  real  parts  gives 


-1  =  '  +  4/3y" 

a  =  tc 1  +  4/3m' 

(2.4) 


These  definitions  of  Q  ^  may  be  related  to  relaxation  phenomena  through  the 


standard  linear  elastic  solid  (Mason,  1958).  The  model  may  be  visualized  as  a 


30 


spring  of  stiffness,  ,  in  series  with  a  parallel  combination  of  dashpot 

(viscosity  n  )  and  spring  (stiffness  M2  ).  At  high  frequencies,  or 
upon  sudden  application  of  stress,  the  system  has  an  instantaneous  elastic 
response  controlled  by  ,  the  unrelaxed  system  elastic  modulus.  At  low 

frequencies,  or  upon  application  of  a  constant  stress,  the  system  modulus 
is  that  of  the  series  springs 

“l  M2 

The  anelastic  dashpot  controls  the  frequency  range  over  which  the  change  in 
modulus  takes  place,  and  therefore  the  frequency  band  over  which  absorption 
takes  place. 

The  stress  strain  relation  for  a  linear  elastic  solid  is  (Liu  et  al., 

1976) 


0  +  T0  O  =  MR(e  e)  (2.5) 

where  *  r)(M^  +  M2)  is  the  stress  relaxation  time  for  constant  strain, 

and  =  n/M2  is  the  strain  relaxation  time  for  constant  stress.  For  a 

•  • 

harmonic  wave,  o  =  juxJ  and  e  «  ju>e  so  that 

o(l  +  jWTa)  =Mr  e(l  +  jure)  (2.6) 

Thus  the  complex  modulus  is 


a 

e 


(2.7) 
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Finally,  Q  is  the  ratio  of  imaginary  part  to  real  part. 


<f 1  -  CO^ 


(2.8) 


where  T  *  /r  x^  is  the  average  system  relaxation  time,  and  where  C(M^  ,  M^) 
is  a  constant  depending  upon  the  elastic  moduli  of  the  system.  In  2.7,  all 
of  the  frequency  dependent  behavior  is  isolated  inside  the  square  brackets. 

That  factor,  which  shall  be  referred  to  as  R^(u>)  is  a  peaked  function  with 
a)  and  (i>  *  behavior  at  low  and  high  frequencies,  respectively. 

These  results  were  extended  to  a  spectrum  of  relaxations  by  Liu  et  al. 
(1976).  They  assumed  a  linear  relation  between  and  x  of  the  form 

1  -  T0/Te  "  C<<1  where  C  is  a  constant.  Then,  the  distribution,  say  of 
,  completely  specifies  the  behavior  of  the  medium.  Liu  et  al.  chose  the 
distribution  function  to  be 


D(Te) 


f1/Te  •  Tl  >  Te  >  T2 
VO  ,  X  >  T  »  X  <  X 


(2.9) 


e  —  1  e  —  2 


where  x^  and  T ^  were  placed  outside  the  frequency  band  of  observation. 
When  these  assumptions  wexe  put  into  the  Boltzman  aftereffect  equation  and 
integrated,  Q  ^  was  found  to  have  the  form 


■i  .  r  -i  u(Ti  *  v  l 

*  C  tan  - p - 

1  +  u  Ti  T2 


(2.10) 


Again,  C  is  a  constant,  and  the  frequency  dependence  has  been  isolated  in  the 
square  brackets.  This  function,  which  shall  be  referred  to  as  l^Cui)  ,  also 
has  w  and  u >  ^  behaviors,  but  about  a  constant  Q  *  between  and 

The  Model. 

The  estimation  of  high-frequency  body-wave  spectra  is  critically  dependent 
upon  the  anelastic  attenuation  model.  In  this  section,  a  modulation  technique 
will  be  developed  which  superimposes  a  frequency  dependence  upon  an  otherwise 
frequency  independent  Q  ^  model. 

From  equations  2.1  and  2.3,  the  amplitude  of  a  propagating  stress  wave 
in  a  viscoelastic  medium  damps  according  to 

dix 

.  -k"x  i(K’x  -  wt)  _  .  -2Qc  i(<’x  -  wt)  (2.11) 

A  6  ©  A  e  6 

o  o 

where  Aq  is  an  initial  value,  x  is  distance  along  the  raypath  and  c  is 
the  phase  velocity.  Since  Q  *  is  a  per  cycle  characterization  of  absorption, 
it  is  convenient  to  integrate  the  attenuation  over  the  raypath  (Carpenter, 

1967) 


=  expf  )  <2-12> 

The  numerical  integration  requires  a  velocity  model  to  obtain  the  travel  time, 

T  ,  and  a  Q  *  model  to  obtain  the  effective  or  average  ,  Q_1 

E 

The  ratio  T/Q^  will  be  noted  as  t  in  the  following  discussion  with 
subscript  a  or  8  for  P  or  S-wave  values,  respectively. 


33 


The  modulation  proposed  here  will  be  applied  to  the  attenuation  exponent 

it 

in  2.6.  Three  simple  assumptions  were  made:  (a)  t  is  a  constant  as  a 

it 

function  of  eplcentral  distance  and  may  be  represented  by  t  =  1.0  and 

■k  k 

t0  =  4.0  for  crustal  earthquakes,  (b)  The  base  t  are  functions  of  depth 

p 

only.  This  is  equivalent  to  assuming  radially  symmetric  elastic  and  anelastic 
Earth  models.  (c)  The  frequency  dependence  of  Q  ,  should  one  be  discovered, 
has  a  multiplicative  separability  from  the  depth  dependence.  That  is, 

t*(z,  f)  =  t*(z)  R2(f>  (2.13) 

where  z  is  hypocentral  depth.  This  last  assumption  preserves  the  simplicity 
of  the  anelastic  attenuation  correction  represented  by  equation  2.12.  Other¬ 
wise  a  separate  correction  would  have  to  be  computed  for  each  frequency  at 
each  distance.  The  modulation,  R2(f)  >  will  be  taken  from  equation  2.10. 

The  modulation,  R^(f)  ,  defined  by  2.8  has  been  tested  but  was  not 
as  successful  as  R2(f)  .  82(f)  models  only  a  single  peak.  If  the  right- 

hand  side  of  the  function  is  used  to  model  the  decay  from  a  spectrum  of  peaks 
the  decay  is  too  fast  for  any  reasonable  distribution  of  relaxation  times. 

The  relative  sizes  of  and  Qg^  may  be  derived  from  equation  2.4 

and  Poisson's  ratio,  a  .  Note  that  both  P  and  S  waves  are  attenuated  by 
shear  loss  mechanisms  (y''  ^  0)  ,  but  P-waves  suffer  additional  attenuation 
from  bulk  losses  whenever  k''  ^  0  .  Bulk  losses  are  generally  neglected 

(Anderson  and  Hart,  1977),  and  setting  k” 


=  0  in  2.4  gives 
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i.3  ,cu2 
,-l  4  'fr 


(2.14) 


For  a  Poisson  condition,  a  =  0.25  and  Q„^/Q  *  =  2.25 

p  a 

*  A 

The  relative  sizes  of  t  and  tQ  depend  again  upon  Poisson's  ratio. 

Cl  P 

If  a  is  assumed  constant  along  the  raypath  and  k' '  =  0  ,  then  P  and  S 

waves  follow  exactly  the  same  path  (Carpenter,  1965)  so  that  QiT^  *  BTg 


fc6  3  a3 


(2.15) 


*  A  i 

For  mantle  velocities,  tQ/t  *  4.16  -*■  4.84 

p  oi 

An  important  feature  of  tg/tct  is  that  this  ratio  of  ratios  cannot 

increase  as  a  function  of  frequency.  If  all  losses  are  in  shear,  then 
*  *  -i  -i 

tg/ta  is  a  constant  because  Qq  ,  Qg  ,  and  Tg  all  change  in  constant 

a 

proportion.  If  a  bulk  loss  operates  over  some  frequency  range,  then  t 

can  only  increase,  because  increases  much  faster  than  the  slight  decrease 

in  due  to  body-wave  dispersion. 

-1  * 

The  absolute  values  of  Q  and  t  are  more  difficult  to  determine. 


The  basic  model  analysis  reported  in  this  paper  used  t  =  1.0  ,  following 

A 

Carpenter  and  Flynn  (1965)  and  tg  =  4.0  (Helmberger,  1973).  The  recent 
Q  *  models  of  Anderson  and  Hart  (1977)  support  a  distance  independent 

A  A 

t„  *  1.0  ,  but  suggest  tfl  >  4.0  with  a  more  distinct  dependence  on  distance, 

a  p 

A  A 

tQ/t  =  4.0  will  be  understood  as  a  minimum  value  for  that  ratio, 
p  a 
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Figure  2.1  shows  the  modulation  of  the  total  anelastic  attenuation  correction 

according  to  the  distribution  proposed  by  Liu  et  al.  (1976)  in  equations  2.9 

and  2.10.  The  curves  marked  "unmodulated"  give  the  correction  defined  by 

* 

equation  2.12  for  a  frequency  independent  t  (i.e.,  frequency  independent 

Q  ^).  The  other  curves  result  from  modulation  by  1^(0))  (equation  2.10) 

where  *  2000  sec,  and  is  as  labelled. 

* 

t  will  depend  upon  hypocentral  depth,  since  earthquakes  may  occur 

below  most  of  the  highly  attenuating  upper  mantle.  Many  authors  (Solomon, 

1971;  Anderson  and  Hart,  1977)  point  out  that  the  anelastic  attenuation  is 

dominated  by  the  asthenosphere .  An  extreme  case  was  suggested  by  Helmberger 

(1973),  who  modelled  a  50  km  thick  asthenosphere  in  which  =  50  ,  surrounded 

by  a  high  Q  (2000)  mantle  and  crust.  Burdick  and  Helmberger  (1974)  suggest 
*  * 

t  =  . 55  and  t„  ■  2.2  for  hypocentral  depths  of  about  600  km.  The  new 

0L  P 

Q  1  models  of  Anderson  and  Hart  (1977)  suggest  a  much  smaller  change, 

amounting  to  only  30%  at  600  km.  For  the  deep  earthquake  reported  in  this 

*  * 

paper,  t  was  set  at  0.6,  and  t„  =  2.4  . 
a  g 

The  anelastic  attenuation  correction  curves  appropriate  for  hypocentral 

depth  of  600  km  are  shown  .a  Figure  2.2,  with  the  modulations  by  equation 

2.10.  The  change  in  the  relaxation  parameter,  T ^  »  will  be  a  subject 

of  the  next  section.  For  now,  note  that  the  reduction  in  predicted  attenuation 

reduces  the  resolution  of  the  modulation  technique  modelled  here. 

The  model  proposed  here  arises  out  of  physical  considerations  of  absorption, 

but  it  also  meets  a  consideration  of  energy  conservation.  S  ..  ffically, 

conservation  of  energy  requires  an  amplitude  spectrum  to  decay  least  as 
-3/2 

fast  as  w  as  u)  ->  oo  .  But  the  exponential  corrections  implied  by  the 

unmodulated  curves  of  Figures  2 . 1  and  2.2  will  raise  the  observed  spectral 


FREQUENCY  (Hz) 

Total  anelastic  attenuation  correction  vs.  frequency  for  crustal  earthquakes.  The  modulated  curves 
are  based  upon  the  model  of  Liu  tt  al.  (1976). 
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Figure  2.2:  Total  anelastic  attenuation  correction  vs.  frequency  for  deep  earthquakes.  The  modu!- 
are  based  upon  the  model  of  Liu  et  al.  (1976). 
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decay  slopes  by  ever  increasing  amounts.  That  is,  for  the  unmodulated  curves 
to  be  appropriate,  the  uncorrected  station  spectra  must  decay  at  a  rate  faster 
than  exponential,  and  such  a  decay  is  not  observed. 

The  modulation  proposed  will  thus  be  interpreted  and  used  as  a  decay 
slope  modification.  A  characteristic  of  each  modulated  curve  is  that  an 
inflection  point  exists  beyond  which  the  curve  is  concave  down  rather  than 
concave  up.  In  the  limit  of  high  frequency,  each  curve  approaches  a  constant 
as  a  function  of  frequency,  suggesting  that  all  frequencies  in  that  range  are 
attenuated  by  the  same  amount.  If  a  station  spectrum  with  an  original  decay 

_ Q 

of  a)  is  to  be  interpreted  according  to  Archambeau*s  (1964)  relaxation 

seismic-source  model,  for  instance,  then  all  of  the  frequencies  in  the  decay 

-2 

slope  must  be  corrected  equally.  To  obtain  an  a)  decay  slope  on  the  same 
hypothetical  spectrum,  a  slope  change  of  must  be  given  by  the  correction 

function. 

Though  the  expected  decay  slope  for  a  given  earthquake-station  pair  is 

-2 

not  known,  the  range  of  slopes  predicted  by  seismic  source  theory  is  u 

_3 

to  oj  .  Each  of  those  slopes  will  be  sought  as  a  function  of  the  relaxation 
parameter,  T2  »  and  the  implications  will  be  examined.  That  is,  a  spectral 
shape  will  be  assumed,  and  Earth  properties  will  be  investigated  under  that 
hypothesis . 

Data . 

Twelve  central  Asian  crustal  earthquakes  and  a  deep  earthquake  in  each 
of  Fiji  and  Brazil  have  been  tested  in  this  study.  Examples  of  the  results 
are  given  in  Figures  2.3-2. 5.  Each  spectrum  presented  is  an  average  of  eight 
or  more  stations  in  the  distance  range  30-60°.  The  individual  spectral  were 
estimated  by  the  autocorrelation  technique  of  Lundquist  (1977).  In  each  case, 
the  uppermost  spectra  were  corrected  by  a  frequency  independent  model  CIT  208, 


LOG  AVERAGE  SOURCE  SPECTRUM  (cm-sec) 


(TIEN  SHAN)  EVENT  OF  (70-06-05) 


LOG  FREQUENCY  (Hz) 


Figure  2.3:  Body  wave  spectra  for  the  event  of  70-06-05.  From  top  to  bottom, 

average  station  spectra  were  corrected  for  anelastic  attenuation  by 
model  CIT  208,  CIT  11  CS2-qM  with  no  modulation,  CIT  11  CS2-QM 
modulated  to  uT ^  decay  rate  and  CIT  11  CS2-QM  modulated  to  u >~3 
decay  rate. 


LOG  AVERAGE  SOURCE  SPECTRUM  (cm-sec) 


(TSAIDAM  BASIN)  EVENT  OF  (63-04-19) 


LOG  FREQUENCY  (Hz) 

Figure  2.4:  Body  wave  spectra  for  the  event  of  63-04-09.  From  top  to  bottom,  the  average 
station  spectra  were  corrected  for  anelastic  attenuation  by  model  CIT  208, 

CIT  11CS2-QM  with  no  modulation,  CIT11CS2— QM  modulated  to  u>~2  decay  rate  and 
CIT11CS2-QM  modulated  to  w 3  decay  rate. 


LOG  AVERAGE  SOURCE  SPECTRUM  (cm -sec) 


HJi 

(FIJI)  EVENT  OF  (72-03-30) 


LOG  FREQUENCY  (Hz) 


Figure  2.5:  Body  wave  spectra  for  the  event  of  72-03-30,  From  top  to  bottom,  the  average  speotr 
were  corrected  for  anelastic  attenuation  by  model  CIT  208,  CIT11CS2-QM  with  no 
modulation,  CIT11CS2-QM  modulated  to  U)~2  decay  rate  and  CIT11CS2-QM  modulated  to 
W3  decay  rate. 
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while  the  second  set  of  spectra  were  corrected  with  model  CIT  11  CS2-QM. 

The  values  of  t  are  given  on  the  figures  and  correspond  to  the  discussion 

of  the  last  section.  The  last  two  sets  of  spectra  in  each  figure  were  selected 

-2  -3 

for  decay  rate  characteristics  of  u)  or  0)  by  specifying  a  frequency 
* 

dependent  t  and  adjusting  .  The  results  presented  are  also 

characteristic  of  the  other  earthquakes  studied. 

Table  2*1  lists  several  important  parameters,  both  of  the  attenuation 
correction  and  the  resulting  interpreted  seismic  source  properties. 

Discussion. 

Several  observations  may  be  made  from  2. 3-2. 5.  First  and  most  obvious 
are  the  problems  with  the  standard  corrections  which  nucleated  this  study. 

Though  CIT  208  provides  an  acceptable  correction  for  the  P-wave  spectra, 
neither  of  the  standard  models  adequately  accounts  for  the  attenuation  of  high 
frequency  S  waves.  The  deviation  from  a  constant  decay  rate  is  not  a  noise 

# 

problem.  Rather,  the  overcorrection  results  from  the  exponential  correction 
function  overpowering  the  decay  rate  of  the  station  spectra.  Unless  the 
station  spectrum  decays  faster  than  an  exponential,  the  overcorrection  must 
occur  at  some  frequency,  whether  noise  is  present  or  not. 

-3/2 

As  mentioned  in  the  last  section,  any  fall-off  less  rapid  than  u>  violates 

conservation  of  energy  in  the  limit  as  w  -*■  00  .  Though  these  figures  demonstrate 

* 

that  an  improvement  may  be  achieved  by  adjusting  the  base  t  values,  the 
improvement  persists  only  over  a  small  frequency  range.  Some  modulation  is 
required  to  reduce  the  exponential  rate  of  Increase  in  the  correction  function 
over  all  positive  frequencies. 

A  very  interesting  result  of  this  experiment  is  the  fact  that  the  P  and 
S  wave  spectra  from  shallow  earthquakes  need  different  modulations  to  achieve 
the  same  decay  slope.  At  first  glance,  one  might  attribute  the  difference  to 
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Table  2.1 


Event 

T/Q 

Slope 

Relaxation 

Seismic  Moment 

Corner 

Stress 

Adjustment 

Parameter 

(10^5  dyne-cm) 

Freq. 

(Hz) 

Drop 

(Bars) 

P 

S 

T2p 

T 

2s 

P  S 

P 

S 

P 

S 

70-06-05 

.42 

1.68 

none 

10.0  3.2 

.16 

00 

o 

• 

6 

1.3 

Tien  Shan 

42.48N 

1.0 

4.0 

none 

00 

o 

• 

eg 

H 

• 

0.7 

4.5 

78.76E 

h  *  20  km 

1.0 

4.0 

“Z 

.08 

.18 

.17 

.24 

7 

36 

1.0 

4.0 

of3 

.16 

.26 

.23 

.25 

17 

41 

63-04-19 

.42 

1.68 

none 

7.0 

13.0 

.20 

.044 

8 

0.9 

Tsaidam 

1.0 

4.0 

none 

.095 

.045 

0.9 

1.0 

35. 7N 

96. 9E 

1.0 

4.0 

—z 

Cl) 

.08 

oo 

t— 4 

.24 

.09 

14 

8 

h  *  33  km 

1.0 

4.0 

CD'3 

.16 

.26 

.28 

.12 

22 

18 

72-03-30 

0.25 

1.0 

none 

69.0  40.0 

.16 

.18 

30 

133 

Fiji 

0.6 

2.4 

none 

.13 

.22 

15 

242 

25.75 

o 

179. 4E 

0.6 

2.4 

— Z 

0) 

.06 

.10 

.16 

.22 

30 

242 

=  532  km 

_  o 

0.6 

2.4 

J 

0) 

.16 

00 

r-i 

• 

.23 

.25 

83 

356 

A  list  of  the  events,  anelastic  attenuation  model  parameters  and 
interpreted  seismic  source  parameters. 


TABLE  1: 


the  different  t 


but  careful  consideration  shows  that  if  both  wave  types 


are  attenuated  by  the  same  mechanisms,  then  the  modulation  must  be  the  same. 

* 

Figures  2.6  and  2.7  show  the  t  (f)  for  the  two  decay  slopes  tested, 
and  both  figures  support  the  same  conclusions  for  shallow  earthquakes:  (a) 
Either  the  high  frequency  cutoff  or  a  window  in  the  absorption  spectrum  may 
be  resolved  for  both  P  and  S  waves,  (b)  The  P  waves  are  attenuated  to  higher 
frequencies  than  are  S  waves.  Since  both  wave  types  suffer  shear  losses  in 
constant  proportion,  the  additional  attenuation  of  P  waves  must  result  from 
a  bulk  loss. 


This  last  possibility  was  discussed  under  the  Theory  section.  The  ratio 

ta/t  must  either  be  constant  as  a  function  of  frequency  or  decrease  relative 
p  a 

to  a  frequency  range  over  which  the  imaginary  bulk  modulus  is  zero.  If 


tD/t  decreases  by  a  factor  of  6  >  0  as  a  function  of  frequency,  then 

p  ex 

equation  2.4  gives 


k' ' 

y" 


I  (6  -  « 


(2.16) 


Under  the  assumptions  of  this  model  and  the  assumption  that  P  and  S  spectra 
should  have  the  same  decay  slope,  the  only  conclusion  is  that  bulk-loss 
mechanism  operates  over  a  narrow  frequency  band  within  the  window  in  the  S- 
wave  absorption  spectrum. 

The  difference  between  modulations  required  for  P  and  S-wave  spectra 
from  deep  earthquakes  is  not  as  great  as  the  same  difference  for  shallow 
earthquakes.  This  could  result  from  either  the  low  resolution,  especially 
of  P-waves,  or  from  the  limited  sample  of  deep  earthquakes  tested.  But  the 
result  may  also  be  interpreted  in  terms  of  Earth  properties.  Suppose  for 
a  moment  that  the  bulk  loss  mechanism  operates  only  in  the  asthenosphere,  while 
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Figure  2  .6:  Observed  t  vs.  frequency.  The  values  of  T£  correspond  to  ^ 
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Figure  2.7:  Observed  t  vs.  frequency.  The  values  of  x2  correspond  to  Figures 
2. 3-2. 5  for  a  high  frequency  spectral  decay  slope  of  uT 3. 


shear  loss  mechanism  operate  throughout  the  mantle.  Then  the  change  in  path 
from  shallow  focus  to  deep  focus  must  affect  the  P  and  S-wave  absorption 
spectra  differently. 

In  Figures  2.6  and  2.7,  the  P-wave  absorption  spectrum  appears  to  be 
reduced  by  a  constant  factor  as  through  approximately  half  as  much  attenuation 
were  seen  by  the  deep  events  as  by  the  shallow.  The  reduction  in  the  effective 
amplitude  of  k' '  would  also  reduce  the  difference  between  P  and  S-wave 
absorption  spectra,  a  result  which  is  also  observed  in  the  figures. 

The  interpretation  of  the  modulation  change  for  S  waves  as  a  function  of 
hypocentral  depth  is  not  clear,  and  further  study  will  be  required  before 
this  writer  will  forward  an  hypothesis.  It  is  likely  that  the  relaxation 
frequencies  operating  in  the  asthenosphere  are  not  exactly  those  operating 
in  the  rest  of  the  mantle.  Though  the  modulation  technique  reported  here,  some 
information  may  be  obtained  in  the  location  of  a  window  or  high  frequency 
cutoff,  but  the  simple  characteristic  of  spectral  decay  does  not  give  enough 
resolution  to  permit  detailed  inversion  for  the  spectrum  of  relaxation 
frequencies. 

A  bulk  loss  mechanism  which  might  operate  over  seismic  frequencies  was 
examined  by  Vaisnys  (1968).  He  argues  that  a  medium  composed  of  a  solid  in 
equilibrium  with  its  own  melt  may  suffer  additional  phase  change  under 
compressive  stress.  If  the  stress  is  harmonic,  the  amount  of  phase  change 
would  vary  with  frequency,  being  peaked  about  to  =  1/x  ,  where  T  is  a 

characteristic  rate  of  the  phase  change.  Essentially,  the  medium  would  tend 
to  freeze  during  one  half  of  the  stress  cycle,  then  thaw  during  the  other 
half  cycle.  Since  more  energy  is  required  to  order  a  lattice  than  to  disorder 
a  lattice,  the  effect  would  be  a  slightly  greater  fraction  of  melt  after  the 


wave  has  passed.  Since  this  mechanism  requires  pressure  changes,  it  would 

attenuate  only  P  waves;  and,  since  the  mechanism  requires  a  partial  melt, 

it  would  operate  primarily  in  the  asthenosphere. 

The  Vaisnys  bulk  loss  mechanism  was  probably  observed  by  Spetzler  and 

Anderson  (1968).  They  observed  attenuation  as  a  function  of  temperature  in 

an  ice-brine  mixture  and  found  a  marked  decrease  in  Qq  for  temperatures 

within  3%  of  melting  temperature.  The  loss  mechanism  is  definitely  related  to 

the  melting  temperature,  T  ,  since  Q  in  the  near  vicinity  of  T 

ta  oi  m 

was  lower  than  Qa  for  either  higher  or  lower  temperatures. 

As  a  side  note,  the  Brune  model  (1970)  stress  drops  were  calculated  for 
the  various  correction  models  tested,  and  the  results  are  given  in  Table 2.1. 

Stress  drop  calculated  that  way  depends  upon  the  cube  of  the  observed  corner 
frequency,  so  small  changes  in  f  can  be  critical  to  the  stress  drop 
calculation.  The  table  itself  speaks  eloquently  of  the  importance  of 

W 

determining  the  proper  anelastic  attenuation  correction.  Until  such  a  criterion 
is  established,  stress  drop  computed  from  body-wave  spectra  must  be  interpreted 
in  terms  of  the  Q  ^  model  used.  In  particular,  comparisons  of  stress  drops 
calculated  under  different  Q  ^  models  are  not  valid. 

Conclusions. 

This  study  has  investigated  the  frequency  dependence  of  Q  *  ,  or  equivalently, 

* 

of  t  .  Evidence  for  such  a  dependence  was  found  by  examining  the  correction 

of  body-wave  spectra  for  anelastic  attenuation.  The  model  reported  here  was 

designed  to  observe  a  high-frequency  cutoff  or  a  window  in  an  otherwise  constant 

absorption  spectrum,  and  such  a  decay  in  attenuation  toward  high  frequencies 

was  observed  for  both  P  and  S  waves.  When  the  criterion  for  the  body  wave 
-3 

spectra  is  an  w  decay,  then  the  spectrum  of  relaxations  controlled  by 
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modulus  has  a  decay  beginning  about  .2  Hz,  while  a  separate  dissipation 

mechanism  operating  on  the  bulk  modulus  maintains  P-wave  attenuation  out  to 

-2 

about  .4  Hz.  When  the  criterion  is  an  ui  dejcay  slope,  then  the  same 
conclusions  apply,  but  the  frequencies  are  about  .35  and  .7  Hz,  respectively. 

The  change  in  the  required  model  parameters  as  a  function  of  depth  de¬ 
finitely  suggest  that  the  bulk  loss  mechanism  is  concentrated  almost  entirely 

in  the  asthenosphere,  while  the  losses  in  shear  are  distributed  in  a  complex 

*  * 

manner  throughout  the  mantle.  Thus  t  and  tQ  must  have  different  behavior 

a  p 

as  a  function  of  depth  over  the  critical  range  of  requencies  of  the  absorption 
window. 

Finally,  the  value  of  stress  drop  computed  from  body-wave  spectra  was 
shown  to  be  very  dependent  upon  the  Q  ^  model  used  to  correct  those  spectra 
for  anelastic  attenuation. 
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i.  SOURCE  THEORY,  STRESS  ESTIMATION  AND  DISCRIMINATION 
C.B.  Archambeau,  Carlos  Salvado,  and  Jeff  Stevens 

Introduction. 

During  the  first  half  of  the  current  contract  period  considerable 
emphasis  has  been  placed  on  the  development  of  very  general  analytical  and 
numerical  models  of  seismic  sources.  To  achieve  the  required  generality 
we  have  considered  combined  analytical  and  numerical  methods.  Specifically 
we  have  considered  Green's  function  representations  of  the  field  obtained  by 
numerical  computation  within  a  relatively  small  zone  surrounding  the  (generally 
nonlinear)  source  of  seismic  energy.  This  allows  us  to  predict  far  field 
radiation  from  three  dimensional  nonlinear  sources  (such  as  earthquakes  or 
complex  explosions) ,  whereas  in  the  past  the  costs  of  numerical  modeling 
could  be  prohibitive. 

A  second  area  of  research  has  been  the  study  of  observed 
earthquakes  using  current  (analytical)  relaxation  models  of  earthquakes.  In 
particular  we  have  studied  events  in  the  Pacific  region  using  m^-Ms  data  to 
estimate  the  tectonic  stress  drops  and  failure  zone  dimensions.  In  this  study 
we  hope  not  only  to  gain  a  better  understanding  of  earthquake  physics,  but 
also  to  verify  predictions  of  m^-Mg  behavior  based  on  relaxation  source  theory 
in  order  to  establish  the  basis  for  event  discrimination.  Further  this 
approach  delineates  those  regions  producing  "anomalous"  (explosion-like) 
earthquakes  and  provides  a  basis  for  the  description  and  understanding  of 
such  events. 

The  following  section  provides  a  brief  description  of  some  of  the  results 
and  conclusions  obtained  from  the  study  of  the  data  for  a  large  set  of 

events  in  the  Pacific  area.  A  final  section  describes,  in  some  detail,  the 
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ntOi-eticai  approach  used  in  the  combined  numerical-analytical  modeling  of 
seismic  sources.  Applications  of  this  theory  are  in  progress  and  will  be 
discussed  in  later  reports.  Current  applied  work  has  progressed  to  the  point 
where  a  specific  source  model  will  be  generated  in  the  near  future  (i.e., 
a  three  dimensional  "ellipsoidal"  earthquake  in  a  non-homogeneous  initial 
stress  field). 

Seismically  Inferred  Nonhydrostatic  Stresses  in  the  Earth’s  Lithosphere. 

The  large  seismic  data  base  consisting  of  measured  values  of  body  and 
surface  wave  magnitudes  from  world-side  earthquakes,  in  the  range  from 
4.5  to  6.5,  has  been  used  to  infer  tectonic  stress  drops  for  events  in  the 
most  active  regions  around  the  Pacific.  The  approach  relies  upon  theoretical 
predictions  of  and  Mg  for  strike-slip,  normal,  and  thrust  earthquakes  at  a 
variety  of  depths,  where  the  rupture  velocity  is  assumed  to  be  a  fixed  fraction 
of  the  local  shear  velocity  and  the  ratio  of  the  rupture  width  to  length  is  also 
assumed  to  be  a  constant.  The  observed  event  data  are  classified  as  to  type  on 
the  basis  of  location  relative  to  known  tectonic  features.  The  stress  drop  and 
failure  zone  dimensions  are  then  determined  from  the  theoretical  predictions  for 
events  of  the  corresponding  type  in  the  appropriate  depth  range.  We  find  that 
the  stress  patterns  inferred  are  highly  variable  spatially,  with  relatively 
small  zones  of  high  stress  near  1  kbar.  The  location  of  such  stress  concentrations 
is  usually  within  a  seismic  "gap"  with  regard  to  earthquakes  much  stronger  than 
those  used  in  the  analysis.  Estimation  of  stress  before  and  after  relatively 
large  earthquakes  shows  that  the  large  events  initiate  in  the  small  zones  of  high 
stress,  with  stress  levels  being  much  lower  in  this  region  after  the  large 
event  but  usually  with  higher  stress  zones  created  near  the  limits  of  the 
failure  zone.  Large  events  nearly  always  show  stress  drops  around 
100  bars,  which  is  much  lower  than  the  level  that  appears  to 


b j  required  toi  initiation.  This  can  be  interpreted  to  mean  that  large 
events  start  within  the  high  stress  zones  but  that  failure  may  be  dynamically 
driven  into  zones  of  much  lower  stress,  sometimes  also  across  these  low 
zones  into  other  zones  of  high  stress.  This  behavior  results  in  a  low 
average  stress  drop  and  the  appearance  of  a  multiple  event. 

Green's  Function  Techniques  for  the  Representation  of  Elastic  Wave  Fields. 

Two  particular  applications  of  the  Green's  function  representation  in 
elastodynamic  theory  to  numerically  generated  wave  fields  are  considered. 

First,  it  is  shown  that  the  Green's  function  representation  of  a  wave  field 
in  the  elastic  zone  outside  of  a  nonlinear  source  region  can  be  used  to 
generate  an  equivalent  elastic  source  that  analytically  reproduces  a  field 
identifical  to  the  numerical  field  in  the  elastic  region.  This  analytical 
representation  of  the  field  can  then  be  easily  employed  to  predict  the 
radiation  at  large  distances  from  the  source.  The  approach  therefore  allows 
truncation  of  the  numerical  computation  after  a  short  time  and  allows  complex 
linear-nonlinear  problems  to  be  attacked  with  much  greater  accuracy  and 
efficiency.  The  second  application  employs  the  Green's  function  representation 
to  effect  a  transparent  boundary  for  the  edge  of  a  numerical  grid.  In  this 
case,  an  infinite  space  (or  half  space)  Green's  function  representation  is 
used  to  predict  the  motion  of  the  grid  points  on  the  boundary  of  the  grid 
based  on  the  motion  of  the  points  on  the  surface  defined  by  the  nearest 
interior  grid  points.  This  approach  causes  the  grid  boundary  to  transmit 
the  field  without  reflection  of  energy  back  into  the  interior  of  the  grid 
system  and  hence  can  greatly  reduce  the  necessary  grid  size  and  improve 
numerical  accuracy.  This  approach  can  greatly  enhance  our  ability  to  obtain 
realistic  representations  of  complex  seismic  sources  wherein  both  nonlinear 
and  linear  phenomena  can  be  modeled  in  an  arbitrary  geometry. 


General  Theory.  We  consider  the  general  linear  elastic 


equation  of  motion* 
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(3.1) 

where  the  stress  tensor  is  given  by 


T=*c  e=C  u  a 

ij  xjki  eki  xjk*  k,£  xjki  u(k,i) 


where  e^  ^  is  the  strain  with  ^  the  symmetric 

part  of  t  a  3U^/3xa  and  is  the  elastic  tensor  with 

symmetry  relations 


cijki  Cjiki  *  Cijik  c  Ckiij 


Both  p  and  the  C^^jj,  may  be  functions  of  the  spatial  coordinates 
but  have  been  assumed  independent  of  time. 


We  may  define  four  vectors  and  four  tensors  which  al¬ 
low  the  equation  of  motion  and  boundary  conditions  for  elas¬ 
tic  problems  to  be  expressed  in  the  compact  general  form:** 


The  general  reference  for  much  of  this  section  is  Archambeau 
and  Minster  [197 V J ,  "Dynamics  in  Prestressed  Media  with  Mov¬ 
ing  Boundaries:  ‘A  Continuum  Theory  of  Failure  (  in  puass, 
Geophys.  J.  R.  A.  S.). 

* 

In  this  discussion,  Cartesian  tensors  were  used  throughout. 
Further,  Greek  letters  used  as  indices  will  always  run  over 
the  range  1,  2,  3,  4  while  Latin  letters  will  have  the 
range  1,  2,  3. 
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L  u  =»  pf 
ay  y  a 


!TTo6neIi 


(3.2) 


a,  8  and  y  take  on  the  values  1,  2,  3,  4  and  the  first  equa¬ 
tion  in  3.2  is  the  equation  of  motion  while  the  second  equa¬ 
tion  expresses  the  boundary  conditions  to  be  satisfied  by 
the  field  0^  on  boundary  surfaces  E.  Here  the  standard 
bracket  notation  for  the  jump  across  a  boundary  is  used, 
i.e.,  for  any  function  F  defined  on  both  sides  E^  and  E^  of 

F  (E  )  -  F(E  )  where  F(E  )  means 
12  1 


a  surface  E ,  then  IF1  £ 

F  evaluated  cn  E  as  approached  from  side  1  and  similarly 


for  F(E  ). 
2 


Here  the  operator  L  in  3 . 2  is  defined  as 

oy 


=  JL 

Jay  "  ax 
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(CafM 


(3.3) 


where  {x_}  is  the  four  vector  (x  ,  x  ,  x  ,  x  ) ,  with  x  the 

8  12  3%  % 

is  the  "elastic-inertia"  ten¬ 


time  coordinate.  Further  C 
sor  defined  by 


ci8Y$ 


Ca8y«  “  “Cijk£;  a'B'Y<S  *  x'2*3 


CoBy«  (Ci4k4  "  C4i4k  *  p6ik?1,k  =  1,2,3 


(34) 


Co6y6  =  0;  otherwise 

where  CagY(S  has  the  same  symmetry  properties  as  did 
that  is: 


Ca8y6  *  ^Bayfi  °  Ca86y  =  Cy5o8 


(3.5) 


In  addition  U  and  f  are  the  "space-like"  four  vectors 
Y  ® 

defined  as 
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{u}  =  (u  ,  u  ,  u  ,  0) 
a  12) 


{f  >  =  (b  ,  b  ,  b  ,  0) 
a  12) 


(3.6) 


while  n0  is  a  four  vector  normal  (in  space-time)  with  a  time 
p  _  * 

like  component  and  is  defined  by  {n0>  =  (n  ,  n  ,  n  ,  w  •  n) 

P  1  2  3  — 

where  w*  is  the  boundary  velocity  vector  relative  to  the 
material  particle  velocity:  w*  =  £  -  V.  Finally,  r^g  is  the 
generalized  inertial-stress  tensor  given  by 


Ta{S  CaBy6  3x^ 


(3.7) 


which  has  the  same  formal  structure  as  the  ordinary  stress 
tensor  T\  ^ ,  but  has  the  added  feature  of  containing  the 
inertial  forces  as  well. 

An  equivalent  form  for  the  equation  of  motion  is, 
from  3.2  and  3.3 


LaB,B 


(3.8) 


so  that  the  equations  of  motion  can  be  expressed  as  the 
divergence  of  the  symmetric  four  tensor  This  allows 

a  Green's  function  representation  for  the  displacement 
field  Uy  to  be  obtained  for  the  general  case  in  which  the 
boundaries  of  the  medium  may  move  with  a  velocity  which 
is  different  than  the  particle  velocity. 

In  the  applications  of  the  Green's  function  solution 
to  be  discussed  here,  we  will  not  be  particularly  concerned 
with  the  cases  in  which  a  boundary  or  boundaries  in  the 
medium  move  with  some  velocity  other  than  the  material 
particle  velocity.  However,  this  feature  of  the  formula¬ 
tion  has  been  used  to  describe  failure  boundary  growth  and 
the  associated  radiation  of  elastic  energy  in  a  prestressed 
medium  by  Archambeau  and  Minster  [1977]. 
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In  any  case,  Archambeau  and  Minster  [1977J  show  that 
Green’s  function  solution  of  the  system  3.2  is  given  by 


4iru^  (x) 


/  Pfa 


GU  d*x 

a  0 


JB  nB  da° 


*  +  f  Id  /  Jw  dv°  I 

1  i  v  4  J 


(3.9) 


where  is  the  concomitant  for  ua  and  the  Green’s  function 
Ga  ~  ^  and  9*ven  by 


J3  Ga  TaB  “  ua  ~aB 


(3.10) 


where 


11  3G 

Gv  =  c  _ 3 

~aB  uaBY5  3x. 


(3.11) 


is  the  inertial-stress  tensor  associated  with  the  Green's 

tensor  G^. 

Y 

From  the  definition  of  <J^  and  the  associated  func¬ 
tions,  we  have  that 


u  f  aG?  3u.  n 

*-  0  0  J 


(3.  12) 


J”  -  0 

4 


where  and  G^  are  Cartesian  tensors  in  the  ordinary  three- 
dimensional  space . 
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Further,  the  various  limits  of  integration  appearing 
in  (9)  are  over  fi,  the  four -volume  in  which  the  time  co¬ 
ordinate  x  5  t  ranges  over  the  interval  (0,  t+)  with 
+  * 

t  =  t  +  e,  with  e  >  0  and  infintesimal  (which  is  used  only 
to  avoid  singular  points  of  the  Green's  function  which  may 
occur) ,  while  the  spatial  coordinates  range  over  the  ordi¬ 
nary  spatial  volume  .  Thus ,  a  v  denotes  the  spatial  sur¬ 
face  of  v  . 

l 

If  we  rewrite  the  Green's  integral  solution  3.9,  using 
the  original  definitions  of  the  various  four  tensors,  we 
can  express  (9)  in  terms  of  the  ordinary  spatial  three  ten¬ 
sors  as 


tT 

4iru  (x,t)  »  /  dt  f  p(x  )  f.  (x  ,t  )G?(x,t;x  ,t  )dv 

m  "  J  o  J  ~ o  k  — o  #k—  Q  o  o 

i 

t+ 

+  fQ  dt.  jf  [GE(*'t;vV  +  «*) 


3G^ 

uk(VV  Jcijkt <So>  «£<*'«=» VV  +<>5t!Swl| 

o  0 


Vv” 


|d[/p(-k‘VV 


3G; 

> 

3t 


dii. 

G™(x,t;x  ,t  )  dv 

K  “  — 0  0  (  o 

0 


(3.13) 


This  result  is  similar  to  the  classical  Green's  function 
solution  except  that  terms  involving  w^,  the  velocity  of  the 
boundary  of  are  present  in  the  integral  over  9V  .  These 
terms  are  negligible  when  the  boundary  moves  with  the  particle 
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velocity,  that  is,  when  the  boundary  is  an  ordinary  material 
boundary,  but  can  be  large  and  comparable  to  the  other  terms 
in  the  integral  over  dv  ^  when  the  boundary  is  a  phase  transi¬ 
tion  boundary.  In  the  applications  of  the  present  work  we 

will  take  all  boundaries  of  V  to  be  material  boundaries  and 

1 

will  neglect  these  terms.  Archambeau  and  Minster  [1976] 
discuss  the  general  problem. 

The  final  term  in  3.9  and  3.13  represents  effects  of 
(generalized)  "initial  values"  of  the  displacement  and 
velocity  fields  and  the  integral  representing  these  effects 
has  been  written  as  a  Stieltjes  integral  since  in  general 
it  has  this  form  when  non-material  boundaries  are  present. 
Thus,  in  this  respect  3.9  and  3.13  are  different  from  the 
classical  representations,  however,  again  it  is  not  neces¬ 
sary  to  use  this  general  form  for  the  applications  here.  . 

For  the  restricted  applications  of  the  present 

discussion,  this  term  reduces  to  the  classical  result  which 

accounts  for  the  ordinary  initial  values  of  and  3uk/3t. 

For  application  of  3*9,  which  has  the  explicit  form 
given  in  3.13,  we  observe  that  the  original  field  equations 
governing  ua  and  are: 


I>  u  —  pf 
ay  y  a 


WM,  =  TaBn8|E 

i  i 


(3.14) 


and 


r 

i 


loth  sets  of  equations  apply  to  the  four-volume  ft.  In  this 
case,  using  the  boundary  conditions  that  apply  for  the 

g 

fields  u  and  G  in  3.9  we  have: 
ci  Y 


4nu 


V 


(x) 


/pf  Gu  d4x 
a  a  o 

ft 


TaBnB 


(3.16) 


This  is  a  rather  simple  form  in  which  all  field  quantities 

on  the  right  side  are  generally  known  since  they  appear  in 

integrals  over  the  medium  boundaries  in  either  space  or 

time  (i.e.,  the  latter  as  initial  values  at  t  =  0  ordinarily). 

When  the  body  force  field  is  negligible  or  entirely  absent 

and  when  the  medium  is  in  equilibrium  initially,  then  only 

the  second  integral  over  the  boundary  of  appears  in 

3.16  and  we  see  from  the  differential  system  in  311#  that 

the  boundary  condition  provides  the  necessary  specification 

of  x  nr\a  on  the  surface  3  V  .  (Here  we  have  used  E  =  3V 
afi  P  1  11 

to  denote  the  boundary  of  v  approached  from  inside  v  while 

is  the  boundary  approached  from  outside  V  ^ . )  In  particu¬ 
lar,  t  0n0  on  3 v  takes  on  the  value  t  0n0jv  which  corres- 

Ct  p  p  1  CL  P  P  I  id 

ponds  to  the  "generalized  traction"  applied  z externally  to 
the  boundary  of  V^.  If  we  neglect  the  boundary  motion  rela¬ 
tive  to  other  effects  as  described  earlier,  then  x  re- 

a6  8 

duces  to  the  negative  of  the  ordinary  tractions,  t^  =  T£kn£ 

on  3V  . 
i 


We  see  from  3.15,  however,  that  we  must  obtain  the 
Green's  function  satisfying  homogeneous  boundary  conditions 
on  the  boundaries  of  ft  in  order  that  3 . 9  reduce  to  the 
simple  form  3.16.  Often  this  is  difficult  and  sometimes,  of 


course,  impossible  for  complicated  problems.  In  this  circum¬ 
stance,-  we  would  want  to  use  an  approximate  procedure  in¬ 
volving  some  Green's  function  that  we  can  obtain  or  is  known 
and  that  is  close  to  that  required  for  a  solution.  Thus, 
for  general  applications  for  which  the  system  3.15  may  not 
be  soluble  in  any  practical  sense  we  would  want  to  consider 
a  Green's  function  given  by  a  system  of  equations  which 
are  "close  to"  3.15,  but  which  are  soluble.  There  are  in 
fact  several  approximate  procedures  that  can  be  followed. 

One  approach  is  to  simply  replace  the  region  ft  over  which 
the  operator  is  defined  by  a  new  region  ft' ,  which  may  not 
contain  all  the  boundaries  of  ft  but  is  larger  in  the  sense 
that  all  space-time  points  of  ft  are  contained  in  ft'  .  Thus, 
we  can  introduce  a  new  region  ft '  which  is  such  that  within 
this  new  problem  space  the  equations  of  3.15  can  be  solved 
for  the  Green's  function.  We  shall  denote  these  Green’s 

Q 

functions  by  T  (x,x  )  in  3.9  where  3.9  is  still  defined 
Y  —  — g 


over  ft  rather  than  ft'  . 
.B 


Since  the  Green's  function 


However,  3.9  will  become 


r*J  are  still  defined  on  ft  then  3.9  has  meaning 
Y  o  p 

when  T  is  used  in  place  of  G  . 

Y  Y 

an  integral  equation  rather  than  a  solution  like  3.16,  and 

g 

while  T*  Will  be  known  and  usually  in  a  relatively  simple  form, 
it  will  now  be  necessary  to  solve  an  integral  equation  for  the 
unknown  field  u^(x>,  by  iteration  methods  for  example. 

Specifically,  we  define  the  class  of  approximating 

O 

Green's  functions  by 


L  rj  (x,x  )  =  (x , x  );  xCft ' 
<*Y  Y - o  a - o  — 


(3.17) 


y“8U,x()  ne  =  0,  s,  & 


where 


<£.•*.)  -  coBy4(x> 


(3.13) 


is  the  inertial-stress  tensor  associated  with  I*v.  We  assume 
,  Y 

that  (2  contains  ft  but  that  at  least  some  of  the  boundaries 

of  ft  are  not  present  in  ft'  (e.g.,  ft'  could  be  an  infinite 

space,  so  that  the  boundary  condition  in  3.17  would  be 

satisfied  at  infinity  and  ft  would  be  contained  in  ft', 

but  would  not  account  for  any  of  the  boundaries  of  ft) . 

O 

Now,  using  3.9  with  T  gives  for  u,  (x)  in  ft: 

•  V  11  — 


4iruvl 


’  /  Ofafa  d'\  -/  dt.  /  ra  [ 
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1  1 


'l  A f< 

t  i 


(3.19) 


where  3V^03V^  is  the  set  theoretic  difference  in  the  spatial 
boundary  of  and  V  ' .  Here  v  ^  is  the  purely  spatial  part 
of  ft  while  v'  is  the  purely  spatial  part  of  ft ' .  Clearly,  the 
third  term  vanishes  if  the  boundaries  of  ft  coincide  with 
those  of  ft.  Here  we  have  taken  the  temporal  parts  of  ft  and 
ft'  to  coincide,  so  that  the  time  interval  for  both  problems 
is  taken  to  be  (0,  a),  or  in  effect  (0,  t+)  due  to  the  causal  properties 
of  the  Green's  functions. 
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We  observe  that  since  ft'  and  ft  do  not  coincide  by 
definition,  then  the  third  term  in  3.19  ,  which  involves  the 
displacement  field  u  on  3  v03v'  ,  is  normally  nonzero.  If 
the  choice  of  ft’  is  a  good  one,  in  terms  of  an  approximation 
to  the  solution  for  u^(x)  using  the  proper  Green's  function, 
then  this  term  in  3.19  will  be  small  or  of  no  importance 

in  the  problem  to  be  solved.  _ _ 

A  good  approach  to  the  solution  of  3.19  when  ft1  has 
been  chosen  to  make  the  third  term  relatively  small  is  to 
take  the  first  iterate  to  the  solution  to  be  the  sum  of  the 
three  other  terms  in  3.19,  which  can  be  computed  directly. 
Thus,  we  take 


(3.20) 


As  the  first  approximation  to  uy  (x) ;  the  second  approximation 
is  then  given  by 


/  "a<V  [^B <*'*,>  nB 

3  V  03  v1 


da 


(3.21) 


with  higher  order  terms  given  similarly. 

These  results  are  essentially  the  only  formal  struc¬ 
ture  that  is  required  to  generate  representations  of  the 
field  which  will  be  useful  in  a  number  of  applications. 

Some  particular  relations  which  amount  to  alternate  expres¬ 
sions  of  the  same  representations  are,  however,  useful. 

The  first  of  these  is  the  expression  of  the  Green’s 
function  representations  for  u  (x)  in  the  frequency  domain. 

By  taking  temporal  Fourier  transforms  of  the  previous  re¬ 
sults,  we  have;  using  3.19  which  becomes  the  exact  solution 
given  in  3.16  when  fl'  **  fi  and  is  the  appropriate  approxi¬ 
mation  for  some  choice  of  H  . 

4*  V<5'“>  "  /  <=Vk  d  \ 

V 

1 

*  I  X  pia 

3V  03V 

1  1 

where  we  have  neglected  the  initial  value  term  involving 

4 

since  in  the  following  applications  we  will  consider  only 
cases  in  which  the  initial  condition  is  an  equilibrium  condi¬ 
tion  and  where  all  boundaries  are  material  boundaries.* 
Further,  since  we  have  taken  all  boundaries  to  be  material 
boundaries,  then  all  tensors  are  three  tensors  and  the  normal 
n^  is  the  ordinary  spatial  normal  on  3v  .  The  notation 
denotes  a  Fourier  transformed  field  quantity  where: 


-  / 

3  V 


rk  tfkt  V  da° 


da 


(3.22) 


* 

Archambeau  [1968]  gives  results  for  the  case  in  which  the 
initial  value  term  is  included,  with  an  application  desig¬ 
nated  to  describe  failure  phenomena. 


i 
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u  (x ,  (u) 
in  — 


u  (X,t)i“XWt  dt 

m  — 


(3.23) 


Here  also  is  the  transform  of  the  two  point  tensor  and 

is  a  function  of  the  spatial  coordinates  x  and  alone, 

while  depending  parametrically  on  oj  (i.e.,  the  source  time 

t^  is  not  present).  Examples  of  these  Green's  functions  in 

transformed  form  will  be  given  in  a  later  section  along  with 

their  time  domain  expression.  Finally,  and  are  the 

negatives  of  the  ordinary  stress  tensors  associated  with  the 

m 


fields  and  respectively. 


In  addition  to  the  direct  expression  of  the  displace¬ 
ment  field  in  terms  of  the  Green's  function  integral  as  in 
3.19  and  3.22,  it  is  often  useful  to  generate  the  displace¬ 
ment  field  from  potential  fields  that  satisfy  simplier  dif¬ 
ferential  equations.  This  can  often  simplify  the  analyti¬ 
cal  problem  since  we  can  use  the  much  simplier  Green's  func¬ 
tions  for  the  potentials  in  an  integral  Green's  function 
relation  that  gives  the  potentials  everywhere  in  terms  of 
their  boundary  values  and  then  obtain  the  displacement  fields, 
stress  fields,  etc.  for  the  entire  region  using  the  dif¬ 
ferential  relation  between  the  displacement  field  and  the 
potentials. 


Specifically,  we  may  define  potentials  xa#  a  **  1»  2 , 
3,  4  such  that  [Archambeau;  1968,  1972]: 


=  v.  !L. .  2v>  [.  !2sl  ♦  f  (3 . 

at2  P  s  pkim  rk  (3*“ 

where  the  four  xa  are  the  three  components  of  the  rotation 
(a  =  1,  2,  3)  and  the  dilatation  (a  =  4)  so  that: 
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y  =  I  e  _ E 

xk  2  ek£m  3xfl 


X  = 
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3U, 

3x, 


(3.25) 


with  the  alternating  tensor.  Equation  3.24  is  just  the 

equation  of  motion  for  the  medium  under  the  conditions  of 
isotropy  and  homogeneity.  That  is,  it  corresponds  to  Eq. 

3.1  when  we  take  the  elastic  tensor,  to  be 


Cijk£  ~  X6ij  6Zk  +  y  (6ii  6jk  +  6ik  5ji) 


(3.26) 


with  A,  y  and  p  in  3.1  to  be  constants.  In  many  applica¬ 
tions  this  is  an  adequate  (local)  representation  for  the 
medium. 

Using  3.24,  it  is  easy  to  show  that  the  individual 
potentials  satisfy  wave  equations  by  taking,  successively 
the  curl  and  divergence  of  this  equation.  We  have,  there¬ 
fore,  that: 
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where 


If  we  consider  temporal  Fourier  transformed  fields,  then  we 
have 
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in  place  of  3.24  and 


Vzx  +  kz  x  **  q 

Ao  o  a 


in  place  of  3.27  ,  where: 


(3.28) 


(3.29) 


Xa  -  ~  /  xa  <x,t)  i-iut  at  (3.30, 

—00 

and  ka  -  u/va  -  (ks,  ks,  ks,  kp). 

Thus,  we  see  that  each  of  the  four  potentials  satis- • 
fied  an  equation  of  the  same  form,  namely  a  scalar  wave 
equation,  and  that  if  we  solve  for  the  xa  then  we  can 
construct  the  acceleration,  velocity,  displacement  and  stress 
fields  from  them.  Naturally,  the  boundary  conditions  on  the 
Xa  are  just  those  given  for  example  in  3.2  and  must,  of 
course,  be  taken  into  account. 

Now  the  Green's  function  integral  representations  for 
the  xa  ate  quite  simple  and  in  addition,  the  scaler  Green’s 
functions,  which  satisfy  differential  equation  of  the  form 

V*G (r,t;r  ,t  )  -  —  —  G(r,t?r  ,t  ) 

-  -0  0-  v2  gt2  -  “0  0 

-  -  4tt<3  (r  -  r  )  6  (t  -  t  )  (3.31) 


is  also  relatively  simple  compared  to  the  tensor  Green's 

function  associated  with  u  . 

rn 

In  particular,  the  Green's  function  integral  giving 
any  of  the  potentials  xo  is  [e»g*r  Morse  and  Feshbach,  1953], 


4irx(r,t) 
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(3.32) 


where  we  have  suppressed  the  index  a  denoting  the  individual 
potentials  since  they  all  have  solutions  of  identical  formal 
structure,  with  v,  G,  x  and  q  having  the  appropriate  index 
depending  on  which  xa  represented,  where  all  initial  value  terms 
for  the  potentials  Xa  have  been  taken  to  be  zero.  Results  applying  when 
initial  values  are  to  be  accounted  for  or  when  moving  boundaries  are 
involved  are  given,  for  example,  by  Minster,  1973. 


Analogously  with  3.22,  we  have  the  representation  in 
the  frequency  domain 


4Trx(r,a))  =  y*  pqG  d v,  +  J  [g  V^x  “  X^Gj  •  n  da0  (3.33) 

v  3V 

vi  i 

where  the  initial  value  terms  are  again  taken  to  be  zero  in 
this  result,  as  was  done  for  3.32.. 

In  both  3.32  and  3.33  ,  we  have  not  introduced  any 
special  notation  for  approximating  scaler  Green's  functions 
corresponding  to  the  earlier  use  of  T™  as  an  approximating 
Green's  function,  since  in  defining  G  and  generating  the 


integral  solutions  we  have  not  made  explicit  use  of  the 
boundary  conditions  for  x  an^  G  on  3  v  .  Hence ,  we  can 
think  of  the  G  appearing  in  these  results  as  being  defined 
on  some  region  ft'Cft,  where  ft’  may  be  identical  with  ft  as  a 
special  case.  In  general,  as  with  the  previous  results  in¬ 
volving  directly,  ft'  will  not  be  the  same  as  ft.  And  the 
Green's  function  in  3.32  and  3.33  will  be  an  approximating 
Green's  function.  Consequently,  the  solution  of  the  result¬ 
ing  integral  equations  for  the  xa  would  preceed  in  the  same 

manner  as  was  indicated  earlier  for  u  ;  where  we  used  an 

m 

iterative  approach. 
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Eigenfunction  Expansions;  Multipole  Field  Representations  of 
Source  Radiation  Fields.  As  a  first  application  of  the  Green's 
function  representations  of  the  previous  section,  consider  the  case 
which  the  dynamic  field  from  a  localized  energy  source  has 
been  computed  by  a  purely  numerical  procedure ,  or  is  known 
by  some  other  means.  For  example ,  suppose  we  compute  the 
motion  of  the  medium  due  to  a  large  explosion,  with  the 
region  in  the  immediate  vicinity  of  the  explosion  behaving 
nonlinear ly  in  response  to  the  high  amplitude  shock  wave. 

At  some  distance  from  the  source  the  stress  wave  will, 
however,  decay  to  the  point  where  the  medium  behaves  linear¬ 
ly.  At  this  point  the  wave  field  will  obey  the  linearized 
equations  of  motion  (1) .  Thus  we  may  view  the  subsequent 
response  of  the  medium  in  the  linear  region  V  ,  exterior  to 
this  nonlinear  zone,  as  being  a  consequence  of  impressed 
field  values  (e.g.,  tractions)  on  the  boundary  of  v  as 
indicated  in  Figure  3.1 (i.e. ,  essentially  this  is  just  the 
Cauchy  stress  principle) .  In  fact,  we  need  not  be  concerned 
rtr  with  what  occurs  in  the  region  V  (Figure  3.1)  but  only  with 
the  values  of  the  field  on  the  boundary  3v  . 

Now  if  we  wish  to  represent  the  field  propagating  into 
due  to  the  impressed  field  on  the  boundary  3Vl»  then  we 
may  use  the  Green's  function  solutions  given  in  Section  I 
directly.  In  particular,  we  can  assume  that  the  initial 

values  of  the  field  in  are  zero,  that  there  are  no  ex¬ 

ternal  time  varying  body  forces,  and  that  we  do  not  have  any 
complications  involving  moving  boundaries  having  rates  dif¬ 
ferent  from  the  particle  velocities.  Then  we  have  for  the 
displacement  field  in  v 
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(T 


AkV 


in 


(3.34) 


Here  T^n^  **  are  the  tractions  and  the  displacements 

on  3  v  ,  which  are  assumed  specified  for  example,  by  the 

numerical  calculation  in  v  and  given  over  the  time  inter- 

2 

val  in  which  they  are  non-zero.  Now  we  may  use  the  infinite 
space  Green's  function  in  3.34  to  calculate  m(x, 

V1  * 

If,  in  particular,,  we  expand  the  Green's  function  r™ 
and  the  traction  derived  from  it  in  terms  of  vector 

spherical  wave  functions,  then  um(k,t)  can  be  expressed  as 
an  expansion  in  terms  of  these  functions.  In  this  case 
um<x,t)  will  be  expressed  in  terms  of  a  vector  multipole 
expansion  about  some  origin  point  which  can  be  taken  arbi¬ 
trarily.  If  we  choose  the  origin  at  a  point  within  v 

2 

corresponding  to  a  natural  symmetry  point  for  the  source  of 

the  radiation  field  we  can  consider  this  resulting  expanded 

wave  field  as  an  equivalent  point  source  which  gives  the 

same  radiation  field  in  v  ^  as  does  the  process  occurring 

in  v  .  (It  does  not,  of  course,  give  the  actual  radiation 
2 

field  within  v  itself,  but  simply  gives  the  analytically 

2 

continued  field  from  v  into  v  .) 

I  2 

The  expansion  of  the  Green's  tensor  in  terms  of 
eigenfunctions  can  be  accomplished  in  one,  two  or  three 
dimensions  for  the  medium  with  boundaries,  so  that  the  region 
v  could  be  layered,  for  example,  and  the  Green's  function 
expansion  would  be  in  terms  of  the  eigenfunctions  for  the 
layered  medium.  Thus  we  can  build  in  the  effects  of  medium 
boundaries  in  the  expansion  for  Um(a,t)  in  v  .  For 
example,  if  the  space  were  a  half  space  or  a  spherical 
region  and  v  (the  "source"  volume)  was  near  or  intersected 
the  outer  boundary  of  the  space,  we  would  use  an  eigen¬ 
function  expansion  for  r1*  involving  half  space  or  spheri¬ 
cal  wave  functions  that  satisfied  the  boundary  condition  on 
the  medium  external  boundary. 
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n 


i-  • 
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Finally  we  note  that  3.34  applies  to  one,  two  or 
three  dimensional  problems  where  3v^  is  to  be  interpreted 
in  general  as  the  boundary  of  the  space,  so  that  if  the 
problem  is  two-dimensional,  then  3v^  is  a  curve.  Then  use 
of  the  Green's  function  appropriate  to  the  dimensionality 
of  the  space  gives  the  representation  of  um# 

We  can  also  represent  um  in  the  frequency  domain 
using  previous  integral  results.  We  have  as  the  frequency 
domain  version  of  3.34  from  3.22  : 


4ttu  (x,w) 
m  — 


dV  ’ 


da* 


(3.35) 


This  representation  can  be  used  in  the  same  way  as  ,  but 
is  somewhat  simpler  to  use  since  the  Green's  function  eigen¬ 
function  expansions  for  the  space  are  easier  to  obtain  and 
express  in  the  frequency  domain  than  are  those  for  the  time 
domain. 

With  regard  to  simplifying  the  analytical  procedures 
required  by  3.34  and  3.34  in  order  to  express  in  v 

as  a  multipolar  (or  eigenfunction)  expansion,  it  is  often 
appropriate  to  generate  the  potentials  x, 
known  values  on  the  boundary  3v 
and  associated  fields  can  be  calculated  in 

i 

Thus,  if  we  choose  to  use  the  xa  for  the  representa¬ 
tion,  then  we  have  that  x„  in  v  is  given  by: 

l 


in  v  from 
Then,  of  course  u 

xn 


4*Xa 


/  «>.  / 

JQ  3  V 


GV  X-XV  G 
o  o 


A 

n 


da' 


(3.36) 


from  Eq.  (32) .  In  the  frequency  domain  we  have 


4*Xa(r,«) 


Gy  X-xV  G 
0  0 


n  da° 


(3.37) 


We  then  can  use  3.27  or  3.28  to  obtain  the  displacement  field. 
Equations  3.36  and  3.37  have  the  same  structure  as  do  the 
more  complex  equations  3.34  and  3.35,  but  are  more  restricted 
in  that  these  results  are  obtained  for  the  case  of  a  homo¬ 
geneous,  isotropic  medium  whereas  3. 34 and  3.35  are  general 
in  this  respect.  Nevertheless,  we  may  consider  the  region 

v  (the  "source"  region)  to  be,  to  a  reasonable  approxima- 
2 

tion  in  most  instances,  confined  to  a  homogeneous  isotropic 
medium  "layer"  or  subregion.  In  this  case  v  is  only  the 
region  within  this  layer  and  3.36  and  3. 37  given  only  the 
field  in  this  subregion.  To  obtain  the  field  in  the  whole 
space  occupied  by  the  (layered)  medium  we  must  adjoin  this 
subregion  field  representation  to  eigenfunction  expansions 
in  the  other  layers  or  subregions;  this  process  of  connect¬ 
ing  subregion  solutions  or  eigenfunction  expansions  is _ 

accomplished  in  the  usual  way  by  applying  boundary  conditions 
at  the  subregion  interfaces  (continuity  conditions) .  Hence 
the  restriction  of  3.36  and  3.37  to  a  subspace  of  the  whole 
space  does  not  generally  reduce  the  usefulness  of  the  expan¬ 
sion  procedure  implied  by  these  results. 

One  important  result  arises  when  v  can  be  considered 
to  occupy  only  one  subregion  of  the  entire  medium  in  which  the 
material  properties  are  essentially  constant  and,  as  a  conse¬ 
quence,  that  ^  can  also  be  considered  to  be  the  extension 
this  homogeneous  zone.  Then  we  can  use  the  infinite  space 
Green’s  function  for  G  in  to  generate  xo»  We  have 

in  this  case,  using  the  Green's  function  expansion  in  spherical 
eigenfunctions  in  a  whole  space: 
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•  i 

Xa(r,«ft)  =  ^  ^  (A£kcosm^  +  B^sinn^J 


1*0  k=0 

•  P*(cos6)  h^2)  (kar) 


(3.38) 


r (2i+i) 

|  4tt 


(COS0 


Here 


C 


Ik 


(2“V 


(1-k) I 
(1+k) I 


We  may,  of  course,  take  3  V  to  be  a  spherical  surface  en- 

2 

closing  v  ,  but  this  is  not  necessary  and  3v  cam  be  any 
2  1 

convenient  surface  enclosing  v  . 

In  case  V2  and  intersect  an  external  boundary  of 

the  medium  (a  free  surface)  then  an  appropriate  appraoch  is 
to  use  a  half  space  Green's  function  representation.  This 
can  be  accomplished  by  using  a  Green's  function  expansion  in 
eigenfunctions  for  the  half  space,  or  by  use  of  an  image 
source,  etc. 

If  we  use  the  time  domain  representation  3.36  rather 
than  3.37  to  represent  the  field,  then 


4*X0(r,t>  -  f  dt  f  I 


G7  x-X^  G 
0  0 


n 


da* 


(3.40) 
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When  the  infinite  space  Green's  function  can  be  approximately 
used,  then  we  have: 


4irXa<£»t) 


6 (t  -t+r*/va) 
8 


V  X 
oAa 


rfi(t  -t+r*/V  ) 

v  T7  0  a 

•  n  i 

xa  V0[  r* 

1 

n 


(3.41) 


However,  we  have,  interchanging  the  integrations  over 

t  and  3V  ,  that 
o  1 


f  1*.  « (t(-t-r*/v0>  V«(W  dtt  = 

i A 


_  i 


r.  V«<Eo't"rVv«> 


/  X“7»( 


6 ( t  -t+r*/v„ )  I 

*-r* - ~J  dt. 


f  a  r«<t,-t+r*/v0)  | 

•J  x0  5pr[ - *-p( - J 


V  r*  dt 

0  9 


£ 


/\  . 

r* 


2 


(r*) 


[-a(trt«*/va)  +ElS'(Vt+£i)  j 


dt. 


r» 

(r*)  ^ 


»t-r*/va) 


t  =t-r*/v 
o  a 


where  r*  is  the  unit  vector  in  the  direction  from  a  source 
point  to  the  receiver  point  (see  Figure  1). 
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Hence,  we  have: 


(3.42) 


Thus,  with  V  x  ,  xn  and  3x„/3t  specified  (numerically  or 
otherwise)  on  3v^  where  this  surface  may  be  chosen  arbi¬ 
trarily  so  long  as  it  encloses  v^,  we  can  compute  x0  in 
in  the  time  domain.  Methods  of  integration  of  3.42 
using  expansions  of  1/r*,  in  terms  of  spherical  eigenfunctions  are 
straightforeward  (e.g.,  see  Archambeau,  1968). 

The  same  procedure  of  integration  over  t  which 
leads  to  3.42  can  also  be  applied  to  3.34,  the  Green's 

v 

integral  representation  of  um(x,t)  ,  using  the  Green's  tensor. 

To  summarize  the  approach  outlined  in  this  application, 
we  observe  that  the  Green's  function  integral  representation 
provides  us  with  a  method  of  representing  the  field  in  the 
elastic  zone  in  analytic  terms  when  the  field  is  specified 
on  the  boundary  of  the  elastic  zone.  This  representation, 
therefore,  provides  us  with  a  way  of  representing  a  non¬ 
linear  source  of  energy  by  a  linear  (elastic)  equivalent 
source  that  gives  us  the  identical  field  in  the  elastic  zone 
surrounding  the  complex  (nonlinear)  source  zone.  Further, 
it  gives  *an  expression  for  this  field  in  an  analytic  form 
which  is  appropriate  for  the  analysis  of  the  further  propaga¬ 
tion  of  this  source  field  in  the  complex  (layered)  elastic 
medium  surrounding  the  source. 
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lication  to  Numerical  Problems:  Transparent  Grid 


Boundaries.  We  can  easily  make  use  of  the  previous  results  to 
show  how  a  numerical  grid  boundary  can  be  made  "transparent" 
to  the  propagation  of  elastic  waves.  We  will  treat  the  simple 
case  involving  potentials  xa  here  to  illustrate  the  method. 

In  particular,  consider  3V^  to  be  the  boundary  de¬ 
fined  by  the  grid  points  at  one  grid  spacing  within  the 

grid  system  as  shown  in  Figure  3.2.  Let  v  be  the  grid  region. 

2 

Now  we  wish  to  treat  the  region  exterior,  to  v  (and  within 

2 

therefore)  to  be  an  infinite  space  so  that  in  this  case 
the  points  on  the  actual  grid  boundary  are  in  v  and  so 
will  transmit  energy  as  if  they  were  in  an  infinite  space 
rather  them  as  the  terminal  points  of  a  numerical  grid.  To 
do  this  we  are  merely  required  to  predict  the  fields  at  the 
grid  boundary  points  B^C  using  the  values  of  the  fields 
specified  on  3v ^  in  an  infinite  space  Green's  function  inte¬ 
gral  representation.  If  we  do  this,  then  the  boundary  of 
the  grid  will  deform  as  if  it  were  in  an  infinite  space 
rather  than  at  the  termination  poihts  of  the  grid  and  no 
energy  will  be  reflected  back  into  the  numerical  grid. 

Using  the  potentials  Xa  (£.»*)  we  observe  that  the 

fields  at  the  grid  boundary  are  given  by  3.42,  where  we 

set  r  =*  r  as  the  coordinate  vector  of  points  on  B  . 

—  — g  g 

Further,  replacing  r*  by  R^  where 

R_  "  |r  -  r  j 
g  — g  — o 


where  r  are  the  "source"  points  on  3v  in  3.42,  then  we 

“0  i 

have  for  the  potentials  on  B  ; 


GRI 
BOUNDARY 

Bn 


oil! 

IHI 

Pi 


Figure  3.2:  Numerical  grid  over  the  region  V ^  with  corresponding  to 

the  surface  at  one  grid  spacing  inside  the  grid  boundary. 
is  viewed  as  the  boundary  to  an  external  infinite  space,  where 
fields  are  defined  by  a  Green's  function  integral. 
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